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ABSTRACT 


NUMERICAL REFRACTION AND 
DIFFRACTION FOR COASTAL AREAS 


by 


ROBERT DALTON SMART 


Submitted to the Department of Civil Engineering on May 19, 1967, in 
partial fulfillment of the requirements for the degrees of Master of 
Science in Civil Engineering and Civil Engineer. 


The purpose of this study was to develop a computer system which could 
be used by engineers in analysis of water wave behavior in coastal areas. 


A wave ray tracing program was developed which uses three-dimensional 
continuous parabolic interpolation for determination of bottom depths and 
derivatives. Tested on analytical plane beaches this method gave results 
which agreed within l per cent of theoretical results. Tested on a 
natural beach, the program agreed favorably with graphical methods and 
numerical methods developed by others, 


For study of diffraction patterns both at ends of semi-infinite 
breakwaters and at breakwater gaps, a program was developed, based on 
Penney-Price methods, which calculates the coefficient of diffraction at 
any point in a diffraction zone. Using formulae derived from the Penney- 
Price methods, a program was developed which calculates directly the di- 
rection of travel of a wave at any point in the diffraction zone, The 
results of this program compare closely with graphical plots of diffraction 
zones which were completed by manual methods. 


These major programs, and a smaller program which locates and traces 
shorelines across an array of water depths, were joined into one system 
through a Problem-Oriented Language. This system reads and works with the 
Same raw data available to an engineer in the field performing analysis. 
It will perform tasks and make calculations as directed by the operator, 
then output all information needed by the engineer in his analysis. 


Either analytical or natural beaches, in feet or fathoms, or combina- 
tions thereof may be input to the program. To simulate tidal fluctuations, 
water level changes may be made at any time during analysis. Breakwaters 
may be superimposed on the array of bottom depths; and, if the breakwater 
or its diffraction zones are intersection by a wave during refraction analy- 
sis, all ray values at the point of intersection are calculated. 


Program output includes, shoreline traces, diffraction coefficients 
with direction of travel of diffracted waves, and wave ray traces. Values 
are calculated at given points along the ray including, among other data, 
direction of wave travel, curvature, coefficient of refraction, and the 
shoaling factor. 
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TERMS AND SYMBOLS 


Wherever possible, terms and symbols used in derivatives and text are 
in accordance with the standards of the Council on Wave Research as pre- 
sented by Wiegel (iglics Hardware and program language requirements prevented 
use of all symbols in actual programs however and the need for additional 
variables may have resulted in notation which is in conflict with the 
standards. Comments within program summaries explain correlation of program 
to text notation. Some symbols have different meaning in different portions 
of the text; however, their use in any particular instance should be clear 


from the text. 


Symbols used in text: 


b Length of wave crest between orthogonals 
C 1) Wave velocity (Phase velocity) 


2) Ina Fresnel Integral, the real portion 


oe Wave group velocity 
D Shoaling Coefficient 
d Depth of water - still level to bottom 
£ Function of one or more variables, e.g. £(X,Y) 
g 1) Acceleration of gravity (32.2 Eee ceta 
2) Function of one or more variables, e.g. go UK.Y) 
H Wave height 
RA Refraction coefficient 


Numbers in brackets refer to References listed in Section VIII. 








Diffraction Coefficient 

20a 

Wave length 

Distance from end of breakwater to point (X, Y) 

1) Increment distance used when tracing a wave ray 
2) Imaginary portion of a Fresnel integral 

Wave period 

Time 

Horizontal coordinate 

Vertical coordinate 

(Delta) change 

(Delta) change 

(eta) surface elevation 

Dies EOD eee tare 

Angular displacement 

Subscript "Oo. refers to deep water value of subscripted 


property 





I. INTRODUCTION 


Waves and their effects constitute a major part of civil engineering. 
Countless books and papers have been written on wave generation, wave 
theory, wave propagation, and wave effects on coastal areas and coastal 
facilities. Waves, and especially ocean waves, are already complex in a 
deep water environment; once the waves enter coastal areas where interac- 
tion takes place with a bottom surface and with coastal facilities such as 
breakwaters, the problems become even more complex. 

In shallow water, wave speed is a function of water depth. As waves 
approach a shoreline, wave crests are bent through the process of refraction 
to conform to bottom contours. This refraction process causes changes in 
both wave height and direction of travel, both of which are essential fac- 
tors in any coastal engineering problem. 

Through the process of diffraction, waves and energy are propagated 
into geometric shadow areas behind breakwaters while reflected waves and 
energy move into regions beyond the ends of breakwaters. Wave patterns in 
diffraction zones are a complex of exponentially decreasing amplitude waves 
spread in an are in the lee of the breakwater, while sets of waves alter- 
nately reinforce and cancel each other in an irregular pattern beyond the 
breakwater tip. 

From the second half of the last century through World War II, these 
problems were solved by empirical methods and through model studies. How- 
ever, due to military necessity for the conduct of operations across beaches, 
extensive research and analysis was conducted during the War, and resulted 
in major advances which formed the basis of the science as it now stands 


today. 








At present, construction of wave refraction diagrams for use in 
engineering analysis requires manual drafting techniques such as those pre- 
sented by Johnson, O'Brien, and Isaacs [2]. However, these methods are 
slow, tedious, and in many ways subjective. Although there have been 
various mathematical solutions for special cases of sloping, circular, 
shelf, and other analytical beaches, these solutions are of limited use to 
engineers in solving problems of natural beaches. 

Penney and Price [3] showed that the diffraction of light is also a 
solution to the water wave diffraction problem. Using Penney and Price's 
theory as a basis, various diffraction diagrams have been published showing 
wave fronts and coefficients of diffraction for standard breakwater prob- 
lems. Due to the extreme complexity of the mathematics of the problem, 
various approximations are usually made to simplify the solutions. Results 
sometimes reflect these simplifications. For engineering studies, diffrac- 
tion analysis is performed through the use of overlays and tracings based 
on dimensionless plots which have been published. These plots are avail- 
able for standard diffraction problems at semi-infinite breakwaters and at 
breakwater gaps of various widths. 

In the interest of good engineering analysis of problems, and in view 
of the capabilities which are now available to the engineer through the 
use of electronic computers, a numerical analysis of these problems which 
is only limited by the basic theory is possible. The basic intent of 
this study was to develop a program system which would accept as input 
the same raw data which is available to an engineer; recognize any limita- 
tions which the engineer wished to impose on the problem; then carry out 
whatever sequence of commands the engineer wishes to use in the thorough 


analysis of any coastal area. 
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The problem of adapting the refraction process to a computer was 
first approached by Lt. G. M. Griswold of the U. S. Navy Weather Research 
Facility [4] in 1963. His efforts and later attempts by others have had 
varying degrees of success. Early programs often gave erroneous results; 
later programs worked better but did not calculate all the information 
needed by engineers; none of the programs were readily suited for use by 
practicing engineers. Analysis of the refraction problem therefore meant 
picking up where prior studies had left off. 

The major problem involved in refraction analysis is that of simula- 
ting bottom topography in a computer. This was the principle difficulty 
experienced with programs which were developed and tested by Harrison [5]. 

In view of the success which had been achieved by Snyder [6] at 
approximating planar curves through a method of continuous parabolic inter- 
polation, it was felt that the same principle, adapted to three-dimensional 
surfaces, offered an excellent yet simple means of surface fitting for 
refraction studies. The method would assure continuous quadratic surfaces 
at any point within the data structure. A method was therefore developed 
for use of three-dimensional continuous parabolic interpolation in tracing 
wave rays over a field of varying depths. The method was tested on slop- 
ing planar analytical beaches and gave results which agreed within 1% of 
the actual theoretical values. When tested on a natural beach, it gave 
results which very closely approximate the graphical procedures now being 
used. Although more extensive testing is desired, the system developed 
appears to offer a means of tracing wave rays which is readily suited to 
use in engineering analysis. The system is of as good or better accuracy 
than present methods, and yet can offer electronic computer capability for 


rapid economical solutions. 
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As best as could be determined, the problem of adapting the diffraction 
process to a computer has not been approached. Various methods have been 
developed for determining coefficients of diffraction using the original 
Penney and Price theories. Without too much difficulty, these methods 
are adaptable to numerical analysis. For accurate engineering analysis 
however, it is also necessary to know the direction of wave travel at any 
given point in a diffraction zone. Although methods are available for cal- 
culating wave phases at any point in a diffraction zone, the only method 
at present to calculate wave direction of travel is through a graphical 
analysis of many points in an area surrounding the point in question. No 
means of directly calculating the direction of wave travel at a point ina 
diffraction zone has been published. 

Solutions of the diffraction problem therefore involved adapting 
existing diffraction methods to a numerical procedure for the computation 
of diffraction coefficients; and the algebraic analysis of a diffracted 
wave front to determine its direction of travel then adapting this analy- 
sis to a numerical procedure. Results obtained from the analysis and 
programs are comparable to those achieved through manual and graphical 
techniques. 

The refraction and diffraction programs which were developed were 
then incorporated along with an algebraic program to locate shorelines in- 
to a system which could readily be used by a practicing engineer. A 
problem-oriented language was used for operation command definition; input 
is the same basic data available to the engineer in the field; output is 


the same product information desired by an engineer in any analysis. 
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Use of this type of program could mean that, given necessary input 
data, engineering analysis of an area could be completed by one engineer 
within a matter of minutes rather than in a matter of months as is the 


case using older methods. 
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Il. REFRACTION ANALYSIS 
A. Refraction Principles 


The velocity of a wave, C, varies with water depth in accordance with 
the basic relationship presented by Eagleson and Dean [7] and others 
21d 


_ gl 4a 
C = ae tanh a (1) 


where g is the acceleration of gravity, L is the wave length and d the 
water depth. With a constant period, velocity and wave length decrease 


as depth decreases. In deep water, 


tanh Sif = 1.00 


and 


Looking at the crest of a long wave moving at an angle to the shore, 
it can be seen that in accordance with equation (1) the deep water portions 
of the crest move faster than the shallow water portions. This velocity 
variation causes the wave length to change and the wave to bend toward 
alignment with the contours. Refraction is the process whereby the di- 
rection of motion changes due to a change in wave velocity. 

Refraction results in changes in wave height, direction of travel, 
and change all characteristics of the wave except for its period which 
still remains constant. The extent of these changes depends on the 


bottom topography, the wave period, and original deep water direction of 


aie 





travel. In performing refraction analysis, it is assumed that as a wave 
travels shoreward, no energy flows laterally along a crest; the energy trans- 
mitted between orthogonals to the wave crests will remain constant. The 
validity of this assumption will be discussed later. Using energy con- 
siderations it can readily be shown, Wiegel [8] (pp. 156), Beach Erosion 


Board [9], and Johnson [10], that 


H = Ho fc 0 | by 
C b 
8 


where H is the wave height. e is the group velocity, and b is the 


orthogonal spacing. The shoaling coefficient, D, is defined by 





C zane 
8 
pe z aes: 4nd/L 
sinh 41d/L 
gees 
n= in deep water 


This shoaling coefficient gives the ratio of the wave height in 
water at some shallower depth to the height of the wave in deep water 


with no refraction effects. 


ase 


This study however will be more concerned with the derivation of the 


coefricient of 2ef naeeron. ir defaned by: 
ON 
b 
eee (2) 
b 


The coefficient of refraction gives the effect of wave crest bending, 
or changes in spacing of adjacent orthogonals, on the ratio of local wave 
height to the deep water wave height. 


The basic principle behind all refraction diagrams is Snell's Law: 


sin Al (Cl 


sin A2 C2 


where A is the angle between a wave front and its respective bottom 
contour and C is as given by (1). This law states that, where bottom con- 
tours are parallel, the sine of the angle between the wave crest and the 
bottom contour is proportional to the velocity of the wave. When all con- 
tours are straight and parallel it makes no difference whether the change 
of depth is a continuous slope or a step, the change in wave length is 
only controlled by the deep water and shallow water wave lengths. 

Although equations have been derived for the solution of circular and 
parabolic beaches by Wiegel [8] (p. 160) and by Pocinki [11], for natural 
beaches, final wave direction is very dependent on intermediate contours 


and incremental procedures must be used for solutions. 


B. Present Methods of Refraction Analysis 


There are two methods in common use today; both are graphical. 
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The Wave Crest Method uses a graduated scale to plot wave advance from 
point to point along a crest. Starting in deep water, successive wave 
crests are plotted until the beach or end of the study area is reached. 

The Crestless Method plots wave orthogonals directly by determining 
shoreward deflections as orthogonals cross successive bottom contours. 

The amount of deflection is obtained from formulae derived from Snell's Law. 

Both of these methods are described in considerable detail in various 
publications including Johnson, et al [2], Bretschneider [12], Dunham [13], 
and Beach Erosion Board [9]. In addition, an apparent improvement to the 
crestless method was developed by Arthur, Munk, and Isaacs [14] which gave 
closer to theoretical results when tested on concentric circular contours. 
Each method has its advantages and disadvantages and the method used in any 
particular situation depends on wave and topographic characteristics, and 
also the engineer/draftsman's capabilities. Both methods, if applied and 
used properly, give results which are in reasonable agreement with aerial 
photographs of study areas as shown by Dunham [13]. 

However, there are two very basic problems with both systems. First 
of all, they depend very much on the operator's skill, ability, and 
judgment. According to Johnson, et al [2], "Since bottom features which 
are comparatively small in respect.:to:.the wave length do not.affect the 
wave appreciably" standard operating procedure is to "use judgement" and 
"smooth out" irregularities. It is left to the engineer to define 
"comparatively small" and "affect...appreciably". 

As shown by Dunham [13] and Pierson [15], different operators can get 
different results using the same system, the same operator can get marked- 


ly differing results using different systems. In particular, the wave 





crest method may smooth out an area of orthogonal convergence and hide 
potentially dangerous situations, (Pierson [15]). 

A second basic problem is that a thorough study of a coastal area is 
extremely slow, time consuming, and thereby inherently expensive. Given 
that very small changes in direction or of period can cause extreme 
changes at a shoreline, it may be necessary to study small increments of 
wave periods for from 2 to 20 or more seconds over approximately 180 degree 
arc ranges, for many areas, at various tide stages. To do a reasonably 
complete job of studying an area can easily require 2 to 3 man months of 
effort; (Pierson, et al [16]). For complex study areas this time can in- 
crease considerably. For smaller low cost facilities or for those in which 
time is of critical importance, thorough engineering analysis quickly 


becomes impossible. 


C. Previous Numerical Methods Used in Refraction Analysis 


In an effort to overcome the limitations and problems of graphical 
wave-ray diagrams, Lieutenant Gale M. Griswald [4], Then of the U. S. Navy 
Weather Research Facility, attempted in 1963 to develop a wave-ray trac- 
ing program based on the numerical methods originally proposed by Munk 
and Arthur [17]. His work, together with that of his associates, resulted 
in programs being written for the CDC 1604, the IBM 7090, and the Bendix 
G-15D computers, (Griswald [4] and [19]). However, these programs gave 
erroneous results in some circumstances apparently because of the bottom 
curve fitting procedures which were used. Results of these programs are 
shown by Harrison [5]. 

This work was then continued by the Army Coastal Engineering 


Research Center. Under contract to CERC, Dr. Wyman Harrison, of the 
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Virginia Institute of Marine Science, developed a method of constructing 
wave rays which appears operational (Harrison [5]). However, for engin- 
eering use this program has many limitations. Because of the linear 
bottom curve fitting techniques used, it is not capable of calculating a 
coefficient of refraction as a ray moves shoreward. Independent passes 
of various programs are required to create input needed by succeeding 
programs, thereby xesulting in rather complex operational procedures. These 
past numerical analyses and attempts to perform refraction analysis on 
computers were milestones. However, they did not go far enough towards 
giving the engineer all the information he needed in a simple efficient 
manner which he would be willing to use in practice. 

The basic approaches and theories behind these programs were sound 
however and optimum use of their methods and experience formed the basis 


mor this study. 


D. Theoretical Considerations for Numerical Refraction 


a) Ray Curvature Equations 





Ray Notation 


Figure l 


= 





As shown in Figure 1, the equations for a ray, S, making an angle, 


A, with the X axis are 


dX dS 


Ae ge Soe (A) 
ts -& sin (A) 
Letting <= = C, we have 
a =C cos (A) (3) 
f= sin (A) (4) 


Munk and Arthur [17] apply Fermat's pene nien to water waves and 
show that the curvature of the ray depends on the gradient of the velocity 


field normal to the ray. Thus 


da_ _ aC 


dt dn 


dA _ dA. dS 


dean csi 


Substituting $2 = C and using 


a 2k ee 
C dn (5) 


© 


da 
ds 


Fermat's Principle states that light travels between two points along 


that path far which the travel time is a minimum. 
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dC 
We can now express the gradient 7 ~ in a more convenient form by using 


the chain rule 


d .3 ax ,a ay 


dn dx dn. dy an 


Noting from Figure 1 that we can show 


cx - sin A 
ey = cos A 
dn 
and 
3 3 3 | 
—_—— soe (Wc of 
= ax (-sin A) 4 ay (cos A) 
Defining a new curvature of the ray, Ros equal to dA/dS we have from 
iy and (6) 
JA 1 dC ac 
SS” = er A —_— 
Aa G50 ee oO ee oe 


It is possible to express wave velocity and its derivatives from 


bottom depth and its derivatives. 


dC _ dd/ax 
dx  ad/ac (8) 
ac _ ad/ay (9) 
dy ad/dac 


One cal calculate dd/93C from equation (1) as follows: 


Ope 
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& 21d 

Cao. tanh (or ) 
Setting 

K' = Ty aa 
kK a 21/gT 

a d 
C = K" tanh 2CK'! 
d = 2CK' tanh ~+ cK") 

= 2CK' x - (2 oe CR) = eek 
dd = i . t! _ 1 _ tt ' ae = = Ke 
tt 

= aca eR? +1n (1+CK") - 1n (1-CK")] (10) 


Given a means of evaluating bottom surface derivatives, equations 
(3) through (10) can be evaluated to find the curvature of a wave ortho- 
gonal at any point and the problem becomes one of tracing this ray of 


Varying curvature from one point to another through a suitable program. 


b) Wave Intensity equations 
Defining a wave intensity factor, or more accurately a ray separation 


maetor, 6, (Beta) by 
B = b/d, (11) 


where b is the distance between adjacent rays, Munk and Arthur [19] de- 


rive the following equation for wave intensity along a refracted ray. 


2 
D 
Dot e eta B= 0 (12) 
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“ - sin A E ee 


c jay) =) 


lI 
| 
QO 
© 
A) 
eae 
Cifr 


p (S) 


2 1 9°c 1 9c 


a x22 SaneAncos 4 le NOY 


Il 
n 
pe 
= 
> 


q(S) 


+ cos’ A (= aC, 
: C 9yY2 (14) 


Factors dd/dC, 9C/¢X, and dC/dY have previsously been defined by equations 
(10), (8), and (9). To evaluate aa O hea. Aa aeaen and 9 C/8Y- we proceed 


as follows. 








ac oc (ody 
a C«~‘d@’:«CY OX’ 
“Cc _ ac . dd) 3d) pod 3*¢ 
9x2 ‘dd2 9x Gx dd 452 


gc _ _ did/dc™ ad? , Bra sax” 

ae 7 ddan os dd/dc 
Similarly 

a*c _ a“ a/ac* (a4)2 , a asay- 

JY2 (dd/dc)> Say AdaG 
and 

9“¢ d-d/dC’ | 8d Bd 9d /9XBY 


IX0Y  (dd/dC)2 9d9Y 3X 9d9d/dC 


To evaluate dady ace we continue the differentiation of equation 


(10). 


ares 


— = K' [1n(1+CK") -—1n (1-CK'"')] +CK' or - Tock"! 








d d a ' Ke z Re Ke Ke! 
acz ~ & EGucKe * ice? tare * xm * 
2 ~ 
CS? + po a)] 
1+CK") Glick) 
2 | 2 2 
d d 7 : Ket Ke a Ge 7 Kat 
acz ~ KR UG Tot + Troe * CCG eRy2 ~ Cecety2?! 9) 


Given a means of computing bottom surface derivatives, at any point, 
equations (13) through (15) can be evaluated. A means of solving (12) 
for 8 is then needed. This is discussed in the succeeding section. 

Once 8 is found, the coefficient of refraction is obtained from 


equations (2) and (11) to give 


(16) 


E. Program Considerations 


a) Basic Program Approach 


Solution of refraction problems involves finding some means of 
tracing a wave ray of known characteristics from deep through shoaling 
water. As explained earlier, at any point on a natural beach, the wave 
characteristics are very dependent on the bottom topography between the 
point and deep water, Procedures which allow for these intermediate 
effects, must be used for solutions. 

Incremental procedures used and tested by Griswald [4] and [19] and 
by Harrison [5] were basically sound. Therefore, similar techniques were 


used in the numerical process developed herein. 
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Using the known wave ray direction ot travel and bottom characteristics 
at a given point, the ray curvature can be calculated using equation (7). 
Using the initial direction of travel, and knowing the curvature, an es- 
timated location of the ray some incremental distance ahead can be 
determined, At this estimated location, it is possible to calculate a 
new ray curvature which can then be averaged with the curvature at the 
initial point. Using this average curvature, one can then revise the es- 
timated new location of the ray. At this newly estimated location, the 
curvature can again be calculated, a new average curvature determined, and 
the process can continue, until the change in curvature between successive 


position estimates reduces to within an acceptable limit. 
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It is possible to calculate all desired wave ray information at this 
new point including the direction of travel. From this new point, the 
next increment of advance can be determined using the same iterative pro- 
cedures. The process can continue until the wave ray runs up onto a 
beach, strikes a breakwater zone, or goes beyond the limits of the array. 

The program developed includes procedures which automatically check 
for shorelines, array limits, and intersection with a breakwater or a 
breakwater diffraction zone. 

To save unnecessary calculations, it is desirable to have the wave 
ray advance rapidly, using long intervals in deep water where bottom 
effects are small, and slowly, using shorter intervals, in shallow water 
where bottom effects are more important. The ray trace process therefore 
includes a scheme whereby, if a long interval is specified for deep water 
calculations, the interval will be reduced in successive stages as the 
ray curvature increases. 

These processes and the manner in which they are used in the program 


are explained in detail in Section IV and Appendix B. 


b) Depth matrix vs. velocity matrix 


Prior numerical methods for making wave-ray computations have been 
based on using an initial computer program to change an original depth 
array into a wave velocity array; then, as necessary, interpolating and 
deriving velocity values and derivatives from this new velocity matrix 
for direct input into the velocity/ordinate equations, such as (7), (13), 
and (14). Making this change to a velocity matrix decreases the amount 
of calculation needed during the processes of tracing a ray and calculating 


refraction coefficients. It does this by eliminating the use of chain-rule 
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multiplication to convert velocity/depth and depth/ordinate derivatives 
into the necessary velocity/ordinate derivatives. 

However, a depth matrix is used in these programs for the following 
reasons: 

(a) An engineer could use mixes of wave periods during the analysis 
of an area without going through the lengthy procedure of resubmitting 
his original depth matrix and converting it to a new velocity matrix which 
corresponds to the changed wave period. 

(b) Water levels can be readily changed to simulate tidal 
fluctuations. 

(c) During possible expansion of the programs for further use in 
coastal analysis beyond refraction/diffraction studies, a depth matrix 


should have fewer limitations than a velocity matrix. 


c) Interpolation surfaces 


Various schemes for representing the bottom topography were tested 
during past analysis of the refraction problem by Harrison and Wilson [5]. 

A "forced-cubic" interpolation scheme was originally used by 
Griswold-Mehr which derived a cubic surface of best fit to the velocity 
values at 12 grid intersections. This surface was forced to pass through 
the 4 values closest to the point in question and was the cubic of best 
fit for the next eight closest eight grid intersections. As a cubic, it 
permitted taking second derivatives for use in wave intensity calcula- 
tions. When this scheme was tested by Harrison and Wilson [5] it was 
found to give completely erroneous results in some cases of a wave ray 


approaching a beach at a shallow angle. 
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A "quadratic-interpolation" surface which was also tested by 
Harrison derived a quadratic surface of best fit to 12 grid intersections. 
This scheme gave good results except when tested within 2 grid units of 
the shoreline. Here it was found to give excessive curvatures. 

A 'linear-interpolation" surface was then developed by Harrison and 
Wilson, and was based on the assumption that velocity values in a grid 
cell could be represented by a plane, This scheme gave good results up 
to within one grid unit of the shoreline but, as a linear surface, did 
not allow calculation of the derivatives which are needed for wave 
intensity calculations. 

For use by engineers, it was felt that calculation of wave intensity 
should be a minimum requirement of any program; this necessitates use of 
at least a quadratic surface. In view of the recent experience of others, 
a quadratic surface of best fit using least squares and possibly with 
heavily weighted central grid values appeared to be a logical method. 
However, this surface would not be continuous as the ray moved from one 
grid to another, Munk and Arthur [17] point out that for successful 
solution of equation (5) a continuous surface is necessary. 

It was found that excellent results at deriving continuous planar 
curves using parabolic interpolation had been achieved by Synder [6] in 
problems not related to refraction. 

Based on the geometrical device of linearly transforming the equa- 
tion of a parabola across an interpolation interval, a smooth curve can 
be developed which is continuous with no abrupt changes in slope even at 
data points. The function and its derivatives at the end of one interval 


are the same as at the start of the next. 
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The method was developed as a compromise between simple interpolation 
with a large input requirement and a more complex interpolation scheme 
with longer computation time. When tested on long segments of a sine 
wave the method gave results which were almost as good as a cubic equation. 
Continuous parabolic interpolation appears to offer a curve fitting scheme 
that is little more difficult than for a simple parabola yet is actually 
a "third degree" equation, Its principal advantage however is the con- 


tinuous nature of the function which results. 
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Continuous Parabolic Interpolation 


Figure 2 


Continuous parabolic interpolation works as shown in Fig. 2. Given 
4 data points, 1, 2, 3 and 4, the problem is one of finding the value at 
abcissa 0 of a smooth continuous curve which passes through all 4 data 
points. In order to calculate the value and slope at abcissa 0, a para- 
bola of the form er bx + c is drawn through points 1, 2, and 3 and its 


projected value at the abcissa of point 4 is determined. An imaginary 
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point, 4°, is located by linearly interpolating the ditterence between the 
projection of parabola 1-2-3 and the actual point 4 value proportional to 
the advance of point 0 from point 2 to point 3. A parabola through the 
imaginary point, 4', is used for calculating the curve value and its 
derivatives at point 0. This interpolated parabola is only valid at 

point 0. It can readily be seen however that as point O moves from point 
2 to point 3 a smooth transition is made from parabola 1-2-3 to parabola 
2-3-4 and a smooth trace of intermediate data values and derivatives is 
generated. 

Although only used in the past for curve fitting, it was felt that 
through the use of multiple parallel and orthogonal parabolas the same 
basic scheme could be adapted to definition of a 3-dimensional surface. 
Just as the scheme generated smooth curves in one plane it should generate 
smooth surfaces in 3 dimensions. 

The simplicity of the interpolating parabola method, especially 
when performed on a computer, and the fact that it met the other refrac- 
tion criteria of providing smooth curves and continuous derivatives at 


all points made it a best choice for a first try at surface fitting. 
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For surface fitting, continuous parabolic interpolation was used as 
shown in Fig. 3. Using the parabolic interpolation techniques just de- 


scribed, 4 independent parabolas were used along vertical planes at A, B, 
2 


Gand D to calculate d, xd and 28 for each of points, (ly YZ ee 
oY 


(3, Y), and (4, Y). A parabola was then interpolated along plane E to 


determine (X, Y) values. One "E" parabola based on 4 depth values gave 


Zz 
a se and 28 at (X, Y). Another interpolated parabola based on the 
od Bis ice oe an interplated ee and ade (X, Y) d imilarl e 
yy Va g ; P ay? AXIY ; and, similarly, on 
based on the 4 2-2 values gave 28 for (Xa). 
oY oY 


Examination reveals that identical results would be obtained by 
reversing the X and Y order and first performing 4 interpolations along 
parallel X planes and the final interpolations along a Y plane. To con- 
firm this fact orthogonal tests were run on real data samples; as expected, 
both tests yielded identical results for (X, Y) and all its derivatives. 

By geometry, it is readily possible to reduce the interpolation 
procedures described above to simple functions. Using the notation of 
Fig. 2, but with the 4 equally spaced data values designated Al, A2, A3, 
and A4; and with D defined as 6/A, which is always measured from A2 


towards A3; one can solve for the value of AQ. 


AQ = AQ & (D/2) (43 — Al) — (0°72) (A = GAG eee 


+ (D?/2) (A& - 343 + 3A2 - AL) 


Its first and second derivatives respectively are 


AO' = (A3 - Al)/2 — (D/2) (A4 - 5A3 + 7A2 = 3A1) 


+ (D*) (AG - 343 + 3A2 - AL) 


AO'' = A3 - 2A2 + Al + (D) (44 - 3A3 + 3A2 = Al) 
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These equations used as functions in a computer program greatly 


facilitate the solution of continuous parabolic interpolation problems. 
d) Solution of the Wave Intensity formula 


Various methods for solving the wave intensity formula equation (12) 


oe Pp PE + dbe= 0 
DS 
are discussed by Munk and Arthur [17]. These methods include: 

(a) Using average constant values for p and q within each interval. 
The equation can be solved and the incremental solutions joined in a con- 
tinuous curve. This solution appeared complex and not readily suited for 
incorporation into the other refraction programs. 

(b) A "WKB" method involves transformation to forms used in quantum 
mechanics, This method appeared more complex than (a) preceeding and 
was ruled out in view of the simpler methods which were available. 

(c) Analogue computer methods could be used to solve this problem 
but were ruled out as beyond the scope of the problem. 

(d) Kelvin's method involves approximating the integral curve by 
fitting together circular arcs. The method is readily adaptable to the 
associated refraction programs. This method was tested by Griswold [4] 
to calculate wave intensity for an analytic field of wave velocity values. 
Wave intensity was then compared’ with actual theoretical values. The 
greatest deviation from theoretical was 1.6 per cent. 

In addition, Griswold [4] used a simple finite difference method. 
Using 3 successive values of 8, designated BN, Bl, and B2, at equal 


Spacing, S, along a curve, one can find by geometry 


ae 





2 
ce = (BN - 2Bl + B2) / ($2) (17) 
dS 


iS 
(ITD 


eu b2e=sBN) /) 2S (18) 


Substitution of (17) and (18) into (12) and solution of the resulting 


equation gives 


_ eee (pom 
BS Gre = Bd Cee 8) BN (19) 


Application of this formula requires that, for three equally spaces 
points, we know Bl, p, and q at the center point, and that we know BN at 
the initial point. The formula calculates B2 at the third point. This 
method requires no iterations, and simply involves projecting one Beta 
value ahead to the next succeeding point, B2. After moving down the ray 
trace, B2 is then used for the next Bl value and Bl becomes BN. This 
method was also tested by Griswold [4] against the same analytic field 
used for the Kelvin method, and was found to give errors of less than 
2.2 per cent. ’ 

In view of the much greater simplicity of the finite difference 
method and the fact that it appeared to give results almost as good as 


this method 
the far more complex, Kelvin's method/was used in the refraction program. 


F. Tests applied to Ray Trace Program 


The theory of the Ray Trace program was tested using two methods. 
First, it was tested on analytical beaches of uniform slope from 
deep water up and onto a beach. These tests could be readily verified by 


means of Snell's Law: 
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sin Al _ Sin A2 
Cl CZ (20) 


and a Beta expression derived from Snell's Law which is given by Munk and 
Arehur [17] 


—_ sin A2 
Sino (21) 


One of these tests, which could be considered a typical case is shown 
in Appendix D. This test shows a wave ray on a plane beach varying in 
depth from 60 feet to zero. The ray has an initial angle of 45 degrees 
and a period of 4 seconds. Its error analysis, based on equations (20) 
and (21) is shown in Table 1. This and similar tests on plane beaches 
gave results which agree within 1 percent of the theoretical values for 


both direction of travel and Beta value. 
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Secondly, the method was tested on a section of natural beach at 
Virginia Beach, Virginia. This test was selected due to the fact that 
previous refraction methods were tested on this beach by Harrison and 
Wilson [5], and results achieved could therefore be tested both against 
graphically constructed wave rays and against other computer programs. 

The most accurate method developed by Harrison and Wilson [5] for 
this beach involved using interpolation with a linear surface of best fit. 
Therefore, the results of the continuous parabolic interpolation tests are 
compared to the linear surface method and the more common graphically con- 
structed rays. The comparison is shown on Fig. 4. 

The accuracy of these comparisons was limited by the fact that the 
input data used by Harrison was not available. It was necessary to es- 
timate input data and his solutions from a 3" by 4" drawing. This could 
have lead to some differences in results but it was felt that data were 
suitable for initial testing. 

These limited results make it appear that the method developed is 
comparable to the graphical scheme and linear surface method. 

It is also important to remember that none of the methods can be 
verified under natural conditions for reasons discussed elsewhere in this 
report. Although graphical methods are most commonly used in practice 
today, different graphical methods would give different results, different 
draftsmen using the same method would give different results. 

One cannot say with certainty which of the three methods most closely 
approximates what would actually happen at this beach. 

In addition, many specialized minor tests were given the refraction 
program to assure that special features such as the breakwater intersec- 


tion check and the gradual reduction of the increment interval do in fact 
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work as programmed. These tests gave what amount to Yes/No results. 
These tests are not enumerated here in view of their lack of significance 


once they work as desired. 


G. Program Limitations 


a. Limitations of the Numerical Methods 


This program is limited in various ways as a result of the numerical 
methods used. Wave ray tracing is based on tracing a wave ray using a 
series of finite intervals. At any point on a natural beach, wave charac- 
teristics are dependent on bottom topography between that point and deep 
water. Any finite interval techniques used may skip over intermittent 
bottom effects. Extensive testing is needed to determine the effects of 
this interval technique. When various intervals were tested on the pro- 
grams developed, the effects were very small; but more testing is needed 
especially on complex bottoms. 

This program is still limited by the ability of parabolic interpola- 
tion to accurately represent the bottom surface, In cases of extremely 
complex bottom topography, it is necessary that grid sizes be small enough 
that bottom features can be reasonably well represented by a quadratic 
surface over a span of 2 grid units, This means that over relatively 
smooth surfaces, large grids are possible; over an irregular botton, 
smaller grids should be used. 

The solution of equation (12) is obtained by a finite difference 
scheme. Although this scheme has been tested against analytical beaches 
by Griswold [4], it should be tested against natural beaches providing 


some suitable means of verifying results can be obtained. 
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b. Theoretical Limitations 


Most of the same original basic assumptions which were involved in the 
graphical construction of refraction diagrams are still present in the nu- 
merical methods. No new theory has been added; the final results are still 
limited by the accuracy of this theory. 

Refraction Theory is based on the assumption that waves are long 
crested, Short crested waves, especially if they cross contours at a high 
angle of incidence, will have refraction coefficients which differ from 
theoretical, Wiegel (8, p 172). 

Although no limits are set, beach slope, roughness, and reflections 
can effect refraction on natural beaches because of their effect on energy 
distribution assumptions used in refraction analysis. 

The assumption that no energy travels laterally along a crest does 
not hold when orthogonals bend sharply or form a caustic envelope. A cer- 
tain amount of energy is transmitted across orthogonals in cases of complex 
bottoms and where sizeable variations in wave height occur along a crest. 
Waves do not actually become infinitely high as indicated by a caustic 
envelope; sometimes they peak and break but usually only a chaotic appear- 
ance results according to Pierson, [15]. Similarly, crossed orthogonals 
may indicate a severence of a wave train and subsequent crossed wave 
crests, but this is an area which is still full of uncertainties. 

Much is still left to the engineer's experience and judgement. 

There is no precise method for smoothing observed contours to that point 
at which they will still represent wave refraction but where neglect of 
minor bottom irregularies is justified. In all cases, accuracy of output 


is limited by accuracy of the depth data available for input. 
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Til. DIFFRACTION 


A. Principles otf Diffraction 


a) General 


Diffraction is the process whereby wave energy is transferred 
laterally along a wave crest and thereby propagated into areas of geo- 
metric shadow behind impervious barriers. It is similar to the diffrac- 
tion of light that takes place when light is transmitted sideways behind 
an object so that the area behind the object is slightly lit. Sound and 
electromagmetic waves experience similar effects. 

Diffraction can best be explained by Huygen's Principle which 
states that each point of a wave front is a source of energy and that 
it radiates outward equally in all foreward directions (Russell and 
MacMillan [18]). The wave motion at any point is the sum of the 
motions induced by all energy sources behind it. In this manner, the 
"end'' of a wave which has been interrupted by a breakwater will act as 
a source of energy and, in the lee of the breakwater, will spread out in 
a circular arc of exponentially decreasing amplitude, (Wiegel [8], 

p. 100). That wave which is reflected off the face of the breakwater 
also has an "end'' which is a source of energy for those regions beyond 
the breakwater tip. This results in two sets of waves, one advancing 
and one radiating, which alternately reinforce and cancel each other so 
as to cause very irregular wave heights in this region beyond the tip. 

As a wave approaches a gap, diffraction takes place in both direc- 
tions as it passes through the gap and two sets of radial waves result 
from the reflected waves. Addition of the advancing waves and two 
radial waves in the lee of the breakwater results in an extremely 


complicated pattern. 





Penney and Price [3] show that the distribution, amplitude, and 
phase of surface waves on a sheet water of any uniform depth can be de- 
scribed by a complex wave function, F(X, Y). Through the use of 
Sommerfeld's solution for the diffraction of light, polarized in a plane 
parallel to the edge of a semi-infinite screen, they are able to resolve 
the water wave diffraction problem. 

The complete solution of the water wave diffraction problem has 
been given in considerable detail by Penney and Price [3], and has since 
been summarized in various forms by many authors, including Putnam and 
Arthur [20], Weigel [21] and [8](p 180), and Bretschneider [12]. The 
reader is referred to these publications for the details and reasoning 


behind equations which are included in summary form only in this study. 


b) Semi-infinite breakwater 


Penney and Price [3] show that the equation for the surface 


elevation, n, in cylindrical coordinates (see Fig. 4) can be stated as 


_ AikC _ikCt 


ies cosh (kd) F(r,8) 


maere A is a constant and k = 27/L. F(+,98) is a function which must 


satisfy 
D, 2 
2 ror 2 Z 
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The diffraction coefficient, K', is defined as the ratio of wave 
height in the area affected by diffraction to the incident wave height 


and is given by the modulus of F(r,@) or 
K' = |F(r,6) | 


os A 





In addition, the wave pattern or phase of the waves can be 
determined from the argument of F(r,8) which is denoted by arg F(r,@). 

The solution of the diffraction problem is therefore the solution 
of equation (22), subject to the boundary condition that the normal 


component of the fluid velocity is zero along the breakwater. 







Xr \ Breakwater 


f 
Ne ee in 


pe Incident Wave Travel 


Notations for Semi~infinite Breakwater 






Figure 4 


Using the notation of Fig. 4, and defining 
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Penney and Price show that in region "S" (i.e., OS O ¢ 80) 
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and in region "Q" (i.e., 8, < 8s (8+) ) 


-ikr cos (8-85) -ikr cos (8-85) 


F(r,80) — e -f(0) e 


+ £(5") eo ikr cos (8+8,) (26) 


where f(co) is defined by 
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uf = + niu? 
f(o) = 5 GL + i) | e 2 du Cy, 


Equation (27) can be evaluated through the use of Fresnel Integrals 


which are defined by 


oe; 
Cea ic -f et" "M2 au 
0 
or alternately, 

CC ) -f cos = TT ie du (28) 
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S(c) -f sin an ac du (29) 
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Substitution into (27) gives 
1 ; 
f(o) = 5 TQ +C+S8S) - i (-C)] 
Using the fact that 
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we can define 


f(-o) = U(o) + iW(o) 


where 
i: 
U = 5 (1 -C - S) 
1 
W = 9 ( S - C) 
setting 
ee eo ikr cos (8-6) £ (0) 
aa a ikr cos (8+6,) £(0") 


f(-o) = Ul +i Wl 


f(-o') = U2 + i W2 


it can be shown algebraically that 
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Ul cos (kr cos (8-8,)) + U2 cos (kr cos (8+6,))) 


i {Wl cos (kr cos (8-85)) + W2 cos (kr cos (8+8,)) 
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Returning to our original equations it is seen from (25) that in 


region 5 


F(r,@) =f +g (39) 


while in region Q, (26) gives 


-ikr cos (8-85) 7 


F(r,8@) =e ft+geg (40) 
In either case, F(r,9) can be expressed as 

F(r,6) =A+ iB (41) 
hence 

Kt = |F(r,6)| = Comey. - (42) 


and the phase difference, $, can be expressed as 
~l 
d = tan (B/A) (43) 


Therefore by evaluating o and o', solving for C and S in (28) 
and (29), and making the transformations indicated by equations (31) 
through (43), we are able to solve for the diffraction coefficient and 


phase difference at any point in the diffraction Zone. 
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c) Breakwater Gap 


Using the following notation 
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Figure 5 


Penney and Price [3] also show that for waves entering a breakwater 


gap at right angles (Gp = 90°) the diffraction pattern is described by 


EO) a= fy see f., af Bo for region S$ (44a) 

F(r,98) = =f tg, + f, + 8, for region T (44b) 

PiGaec) = me ie: -f, +g, for region Q (44c) 
1” 817°" 2 ee 


where the subscript indicates the breakwater wing from which the func- 
tion is measured. 
This solution is only valid, however, for waves entering normal 


to the gap. Upon entering a gap at an oblique angle, a given wave 


b= 








reaches one tip first and diffracts there before the remainder of the 

wave reaches the other tip where it will also diffract. Therefore the two 
wave patterns are out of phase unless the wave crest happens to reach 

the second tip an exact multiple of wave periods after reaching the 

first. Equations (44) do not, in general, give valid results with 

oblique waves. However, an approximate solution for oblique waves is 


possible through the use of techniques which will be discussed later. 


B. Present Methods of Diffraction Analysis 


It is apparent that the solution of the diffraction equations is 
a lengthy procedure when done by manual methods. Fortunately, however, 
they can be solved independently of wave speed and wave period; and 
solved in terms of wave length, L. 

When plotted in terms of X/L and Y/L coordinates, diffraction 
charts become dimensionless plots of diffraction coefficients which 
can then be applied to any water wave diffraction problem in water of 
uniform depth regardless of wave period or wave speed. These solutions 
have been completed and the results have been included with virtually 
every diffraction publication, including Wiegel [8] (p. 180), 
Bretschneider [12], and Beach Erosion Board [9]. For the engineer in 
the field, it then becomes simply a matter of picking one of these 
dimensionless plots, expanding or reducing the scale so as to fit his 
particular situation, and with the use of a properly scaled overlay, 
analyze the particular diffraction situation at hand. Different over- 
lays can be made to represent each wave period to be studied. 

The same procedure applies for studying diffraction at breakwater 
gaps. In these cases, the ratio of gap width, b, to wave length con- 


stitutes a second independent parameter; hence, a series of dimensionless 





plots are required. These have been published by many authors, including 
Johnson [24], and Wiegel (3, p 189), for differing b/L ratios from 1 to 
5. Beyond 5 wave lengths it is felt that the 2 breakwater tips which 
form the gap are far enough apart that they can be considered independent 
for all practical purposes. 

Plots of breakwater gaps for waves approaching at an angle are 
available (Weigel (8, p 193) and Johnson [10]) based on the theories of 
Morse and Rubenstein, which were developed by Carr and Stelzride [25]. 
The calculations for these plots are based on series of Mathieu functions 
which give exact theoretical solutions for incident waves of any angle 
of approach. 

In addition, Johnson [24] showed that a modified gap width as shown 
in Fig. 6 could be substituted for the actual gap and results achieved 
which compare very closely to the exact theory of Morse and Rubenstein. 
This modified gap solution allows the use of Penney and Price equations 


which were developed earlier. 
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Figure 6 
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In areas sufficiently removed from the breakwater, simplifying 
approximations can be made for some terms of the various formulae. 

These approximations have been verified experimentally by Putnam and 
Arthur [20]. However, in regions close to a breakwater a complete 
solution is needed, 

In addition to Putnam and Arthur [20], diffraction theory has been 
verified experimentally by Blue and Johnson [26]. Their experiments 
showed that actual amplitudes agreed very closely with theory in the lee 
of breakwaters and were somewhat less, but within 10 percent of, theore- 
tical beyond the breakwater tip in the unsheltered zone. This results in 
generally conservative results, i.e., overdesigned facilities, if the 
theoretical results are used. 

In addition to coefficients of diffraction, it is also necessary 
to know the direction of wave travel at any given point in a diffraction 
zone. From equation (43), it is readily possible to calculate the phase 
difference between the diffracted wave and the wave front of the original 
incident wave. With this as a basis, it is possible to calculate a 
number of phase differences within a zone, then by connecting zones of 
equal phases determine actual wave crest positions and direction of 
travel. However, no published method has been found which would di- 


rectly give a direction of wave travel for a given point. 


C. Previous Numerical Methods 


No publications were found concerning numerical analysis of water 
wave diffraction problems. However, the present formulation of the 
diffraction problem is readily susceptible to calculation of coeffi- 


cients of diffraction in their present form, hence it is possible that 


0 





portions of the problem have been done by computer as an aid to graphical 


analysis but were not considered worthy of publication. 


D. Theoretical Considerations for Numerical Diffraction 
a) Coefficient of Diffraction, K' 


Except for evaluation of Fresnel Integrals, equations (23) through 
(42) can be solved on a computer in essentially their present form. In 
the programs developed, a completely random orientation of breakwaters 
and incident waves was desired. This and a desire for reasonable program 
efficiency and memory size cause substantial bookkeeping problems but 
create no need for any new methods or procedures. 

As shown in Appendix A, a series solution was developed for the 
evaluation of the Fresnal Integrals, This was used for sigma values 
less than 3.0. For larger values, sine/cosine approximations given by 
Sparrow [23] were used. If used with sigma values greater than approxi- 
mately 3.0, the developed series solution requires greater precision 
than normally used on most computers. Double precision could be used 
but tests indicated that satisfactory evaluation, good to within 0.1 
percent of the Integral, could be achieved by using the two methods. 

A third approximation formula is given by Jahnke and Ende [22] which 
is suitable for evaluation of values in the intermediate range between 
the series solution and the sine/cosine approximations but this was 


not felt necessary. 


b) Direction of Travel of Wave Crest 


As noted earlier in equation (43), the phase, $9, of the diffracted 


wave is given by ¢ = nie (B/A) 
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The equation of the wave crest is given by 
d (z,x) = Const 


The normal to the wave crest is a vector which is given by 
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in which 2 and k were the unit vectors in the x and z directions 


respectively. Therefore, representing the slope of the normal as 


ah we have 
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The partial derivatives in (44) can be expressed in terms of A and 


B. Using (43) this gives 
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Equation 47 is our desired solution for the direction of wave 


travel. 
Given A and B from equations (39) or (40), we must now evaluate 
dA/9X, JB/OX, JA/dZ and 3B/dZ. For these we will need the derivatives 


dr/dX, dxr/dZ, JB/dX, and 96/dZ. For the breakwater along the +Z axis, 
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For the breakwater along the -z axis 
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If we consider X = f£(r,86), the chain rule gives 
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Similarly, one gets 
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Using the same techniques we can evaluate Se and se, but in this case, 


results differ with the breakwater orientation. If the breakwater is 


along the +z axis, 
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while if it is along the -z axis 
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Still using the notation of Fig. 7 the following convenient notations 


are used: 


Keen) Cis 
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Introducing this notation into equations (39) and (41), in region S, one 


has 
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Taking derivatives 
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Similar equations for 
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(48) 


dFP 
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A, B, and their derivatives apply in region 


and may be derived from equations (40) and 
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These equations are not detailed 





in this report in view of the fact that they offer nothing new in the 
way of knowledge and are simply exercises in bookkeeping. 

For evaluation of these last groups of derivatives, one must deter- 
mine Ul, U2, Wl, W2 and their derivatives. Consider first the evaluation 


of Ul and Wl. We recall from (31), (32), and (35) 


ui = - (1-S-c) 
Wi = > (S-C) 


Taking derivatives 
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To evaluate derivatives of S and C it is convenient to set 
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Equation (29) therefore gives 
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Using the same techniques we can evaluate 8 derivatives of S and C and 


get 
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From (23) we can derive 


1 
= Bui ae 
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We evaluate U2, W2 and their derivatives similarly, using é' and its 


derivatives, i.e. 


ao" | E 1 
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Therefore, given values of r, 6, 84> and k, we can back substitute 
through the preceeding series of equations to ultimately evaluate the 


direction of wave travel as given by equation (47). 


c) Alternate Methods Considered 


The two methods used in diffraction analysis by engineers have been 
the Penney-Price and Morse Rubenstein methods. When evaluated by Carr 
and Stelzride [25], both were found to have advantages and disadvan- 
tages but the Penney-Price method was found to be the better. Ever 


since, the Penney-Price has been used almost exclusively except in the 
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evaluation of breakwater gaps where the more exact theoretical solutions 
of Morse Rubenstein have often been used. Principally due to the fact 
that more testing had been accomplished and much more information was 
available in publications on the Penney-Price methods, this approach was 
used in the analysis. 

Once committed to this method, when it was later decided to add 
solutions of the breakwater gap problem to the program package, the 
Penney-Price solutions using the imaginary breakwater gap developed by 
Johnson were used in the interests of simplicity. Although not an exact 
solution when applied in the vicinity of the breakwater wings, in the 
hands of an experienced engineer these areas should cause no real problems. 

Many simplifications which omit negligible terms in various re- 
gions of a diffraction zone could be made. Due tothe unknown effect 
which these approximations would have in the determination of the direc- 
tion of wave travel, an exact solution was used throughout. The actual 
computer time saved by making approximations would be very small and was 
not felt to be worth the chance of having to rewrite programs from 
approximate to more exact solutions if the approximations did not work. 
Working from an exact solution however, approximations could be added 
at a later time and results readily verified against the original 


solution. 


E. Tests applied 


To test the accuracy of these programs, they were used to calcu- 
late diffraction coefficients and directions of wave travel at many 
random points using various breakwater coordinates and various direc- 
tions of travel for the incident wave. Results were then checked 


against graphical solutions which have been published. 
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Problems were encountered in this verification however, due to the 
fact that virtually every published refraction diagram was based on 
approximations of one kind or another (usually not mentioned in the 
write up) which had a substantial effect on the quality of the plot. In 
fact, no two published diffraction diagrams from different authors 
appeared to agreee with each other as a result of these various approxi- 
mations. The program results did not completely agree with any published 
diffraction diagram. However, in the calculation of the value for a 
coefficient of diffraction, good correlation was found between the program 
and the diffraction chart in Wiegel (8, p.183). Also, in the calculation 
of the direction of wave front travel, good correlation was found with 
the diffraction chart published by the Beach Erosion Board (9, p. 38). 
Table 2 compares program results with diffraction coefficients inter- 
polated from Wiegel and with directions of wave travel interpolated 
from the Beach Erosion Board. It should be remembered that in both 
cases values interpolated off the charts are no better than the accuracy 
with which the chart was drawn or one's ability to interpolate off 
them. In addition, analysis of the basic equations, especially in the 
region around and beyond the breakwater tip, where the most discrepan- 
cies are found between published charts, indicates that the results 
achieved are in fact reasonable and in good agreement with the complete 


theoretical solution. 


F. Limitations of Numerical Diffraction Programs 


In that these programs are built on the same theoretical basis as 
the past plots, they suffer the same limitations. No new theories or 
improved results can be claimed over a completely accurate plot. The 
program's chief asset is its ability to produce results in a small 


fraction of the time required by graphical methods. 
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Its output is still no better than its input data. Study of 
complete wave spectra are still needed for thorough analysis of an area. 
The program should however make this analysis easier than it has been 
in the past. 

There are still differences between theoretical and experimental 
results as shown by Putnam and Arthur [20], Blue and Johnson [26], and 
Carr and Stelzride [25], which were mentioned in Section III B. Although 
the agreement between theory and experiment is satisfactory, theory 
does overestimate wave heights by up to 10 percent in the zone beyond a 
breakwater tip. This is probably due to a significant loss in wave 
energy during the process of reflections, (Terry [27]). All calculations 
are based on ideal breakwaters with 100 percent reflection. In practice 
incident waves are only partly reflected and are partly destroyed by 
turbulence (Penney and Price [3]). 

A basic assumption is that the wave height is small compared to its 
length and that the wave profile is nearly sinusoidal. Although the ex- 
periments mentioned above indicated that the results are not seriously 
altered by waves of moderate height, design for solitary waves likely 
to be generated during severe storm conditions could be susceptible to 
error, 

Gap solutions have only been experimentally verified by Blue and 
Johnson [26] with a wave angle of approach less than about 60 degrees; 
however, they are the only solutions available so they are actually used 
in practice over the entire 90 degree range. 

For short breakwaters of the ribbon type, with diffraction about 
both ends, the program solutions would not be valid in a region of dif- 
fraction effects from both ends. In practice however, the breakwaters 
are usually of sufficient length to achieve good results by diffracting 


each end independently. 
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IV. PROGRAM SUMMARY AND ENGINEERING FEATURES 


The program developed was written so as to facilitate its use by 
engineers in the analysis of coastal areas. The following summaries 
give highlights of their operation, Details of the main program and 


all subroutines are explained in the appendix. 


(1) Main Program 


The purpose of the main program is to read data and commands, set 
variables, and call on the various subroutines to take appropriate action 
when necessary. The main program is written so as to read problem- 
oriented language input. In this way, the operator uses the same langu- 
age in describing input data and tasks to be performed as he uses in 
normal every day work. When appropriate, a written output of data is 
provided to assure a complete record of the calculations performed. 

An attempt was made to provide most of the commands that an engineer 
would have occasion to use in an analysis of a coastal area. Although 
this program was, of necessity, written for batch processing, it is well 
suited to console input and random interrogation by an operator follow- 
ing the progress of a solution. 

Since much of the refraction/diffraction analysis process is still 
left to engineering judgment, the flexibility was provided for the en- 
gineer to use his judgment whenever it is deemed necessary. A minimum 
of information is set by the program, it depends on an engineer for 


important variables. 
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(2) Refraction Subroutine 


This subroutine traces wave rays across an array of water depths 
until the ray goes onto a beach, goes off an edge of the array, or 
enters a breakwater zone. During this trace, all information concern- 
ing the wave ray which could be of use to an engineer is output, in- 
cluding location, water depth, speed, direction of travel, shoaling 
coefficient, and coefficient of refraction. 

Due to the fact that in many cases refraction analysis is made in 
breakwater areas, this program stops any wave traces upon crossing a 
breakwater or upon entering a breakwater diffraction zone. At these 
points of intersection, all necessary wave ray information is determined 


and output for use by the operator in further analysis. 


(3) Diffraction Subroutine 


This subroutine calculates the coefficient of diffraction and 
direction of wave front travel in a diffraction zone, This may be a 
diffraction zone at the end of a semi-infinite breakwater, or a diffrac- 
tion zone at a gap between two breakwaters. Given the location of one 
or two breakwaters and necessary data on the incident wave, this progran 
analyzes the diffraction pattern at any specified point and outputs the 


Soetticient of diffraction and the wave front direction of travel, 


(4) Shoreline Tracing 


This subroutine will find and trace shorelines across the array of 
water depths. Although this feature is most useful if the program is 
converted to some sort of a graphical output device, it can also be 


helpful when working in an area where tidal fluctuations can have a 
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noticeable effect on the shoreline location. Once a depth array is 
input to the program, presumably based on a "mean-seqa level", tidal 
changes can be simulated through the input of water level changes; if 
these changes affect the shoreline, this program will determine its new 


location, 


(5) General Features 


To facilitate maximum adoption of this system by any interested 
users, it was written in the FORTRAN IV programming language. Although 
written for the IBM 360 computer, it can readily be adapted to most other 
hardware simply by changing the alphameric character constants in MAIN 
and by adopting the read/write statements to available devices. 

The program is now written for any water depth array size up to 50 
by 40 grid units. This can be changed up or down as desired, simply by 
changing the DIMENSION cards. No other changes are needed. 

The entire system at present requires about 10,000 words of memory 
on an IBM 360 - Model 40 computer. However for smaller cores, it is 
possible to use the program in portions. For example, during refraction 
calculations, subroutines, RSSHOR which locates the shoreline, and 
RSDIFF which performs diffraction calculations, are not needed and mini- 
mum sized dummy subroutines may be inserted for them. During diffrac- 
tion calculations, only RSDIFF and MAIN are needed. For tracing a shore- 
line, only MAIN and RSSHOR are needed, Although this method of operation 
would cut down on some of the operator flexibility, it can facilitate 
program use on hardware too small to accept the entire package and it 
could enable an operator to trade off unnecessary program for an in- 


creased array size. 
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V. POSSIBLE SYSTEM EXPANSION AND IMPROVEMENTS 


The programs developed show some of the possible applications of a 
computer by an engineer in the study of .cvoastal areas. While writing 
these programs it became apparent that many additions and improvements 
could be made to them which would still further facilitate engineering 
analysis. The following are some of these possible expansions. These 
explanations are not intended to be complete analyses but are only to 


point out possible areas of further study. 


(1) Wave Front Tracing 


In many studies, plots of wave fronts can be very useful in 
graphically depicting a situation. The present program, by tracing wave 
rays, will not give any wave front data other than the fact that the. 
wave fronts are orthogonal to the ray and simulated fronts could possi- 
bly be added to the rays by an experienced draftsman or engineer. It 
would be possible however to expand the present program to output wave 
front coordinates. The present ray trace program is based on increments 
of equal distance between points. This could be changed to increments 
of equal time. As the iteration process, which is needed to locate 
succeeding points, take place, an average speed would be calculated and 
the interval between points could be based on this average wave speed 
and the given periods of equal time. If all wave rays were started 
along a wave crest, (this ability to start multiple rays could also be 
added to the program) each successive coordinate point output by each 
ray would be a successive crest position. This feature would also have 
the beneficial effect of reducing the increment interval as the ray moves 


into shallower depths. 
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(2) Graphical output 


The refraction/diffraction problem is ideally suited for a 
graphical output. The combination of graphical output and console input 
by an engineer would be an extremely powerful tool in the analysis of 
any coastal area. Early versions of the program were run on an IBM 
1620 connected to a Gerber Plotter. Unfortunately, the program was 
not perfected prior to removal of this facility from M.I.T. 

Although not used by the author, perhaps the most interesting and 
promising output would be accomplished on an oscilloscope. Through a 
combination of console and light pen input to the program, an engineer 
could rapidly analyse diffraction and refraction problems. With a mini- 
mum of effort he could investigate far more possible solutions than 


would be practical by any other method. 


© Breakers and Surf 


At present, most predictions of breaker and surf locations along 
a shoaling beach are done empirically. These locations can be very 
important in the analysis of a beach due to the fact that refraction 
does not take place beyond the breaker zone and direction of travel of 
the breaker is very important in determining littoral currents. If 
breaking or surf occurs on a reef at some distance from a shoreline, 
especially at different tidal elevations, this can be very important. 
For some studies therefore, it may be desirable to add to the program a 
means of predicting the breaker location. This could be done by 
determining a desired breaking criteria and adding it as another check 


item to the refraction calculations. 
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(4) Refraction by Currents 


In some areas, especially at estuaries, one may have to allow for 
refraction caused by currents. In the past, using graphical methods, 
it was impractical to account for local currents in refraction studies 
and they were either neglected (Dunham [13]) or accounted for by engin- 
eering judgement. A theory for the refraction effects of currents is 
well established and has been published in some detail by Johnson [10]. 
Due to the possible occurrence of currents in coastal areas resulting 
from long shore currents, currents caused by topography, flow from re- 
gions of high refraction to areas of low refraction, eddy currents, and 
rip currents, it may be advantageous to have an ability to account for 
their effects. This would be a desirable addition to the present 


program. 


(5) Output Summaries 


A suitable means of summarizing output is needed. Various methods 
of portraying refraction for a location, showing the effects of period, 
and direction of wave travel have been used in the past by Johnson [10], 
and the Beach Erosion Board [9]. As the volume of data increases how- 
ever, most of these methods become confusing. It would appear that this 
is an area requiring considerable common sense on the part of the engin- 
eer to assure that a detailed study with noteworthy results is not lost 
or obscured by a poor presentation. It would be most desirable, of 
course, to have a suitable data presentation come directly from the 


computer. 
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(6) Approximations for Diffraction Process 


As noted earlier, diffraction theory, when applied to water waves, 
generally gives conservative results. It was also found that in many 
instances, approximations of the complete theoretical solution brought 
results into closer agreement with the tests than did the complete 
solutions, (Blue and Johnson [26]). The Beach Erosion Board [9] points 
out that becausé of non-uniform validity of various orders of wave 
theories, the best approximation is not necessarily the highest order. 
LaCombe [30] proposed approximate solutions to the diffraction problem 
which may offer some avenues of approach. It may therefore be desirable 
to adjust the present program to allow for various approximations. 
These approximations can be compared with the complete solution and can 


be checked against experimental data for an empirical comparison. 


(7) Command structure 


To the maximum extent possible, the program was designed to 
facilitate expansion and improvement. Additional commands may be added 
to MAIN simply by adding their appropriate alphameric characters to the 
listing and adding a block of appropriate instructions to the program. 
The existing commands could be improved through the addition of allow- 
able abbreviations and shorter names. Input is well suited to free 
format and this would be a most desirable feature, especially if console 


input is to be used. 


(8) Multiple Waves 


In many instances, different period waves and waves travelling 


from different directions will arrive at a given point at the same 
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time. There are methods for determining the resultant of multiple 
waves as shown by Eagleson, et al [7]. This is not a simple addition 
problem and, in view of its complexity and common occurrence, it would 


be a desirable feature which a design engineer could use. 


¢ 


(9) Combined Refraction and Diffraction 


All diffraction theory is based on waves of constant speed. In 
practice, breakwaters are usually built in areas where shoaling is an 
important part of the problem and the wave speed variation is important. 
To date, no combined theory of refraction and diffraction has been 
developed. According to Bretschneider [12], and many other authors, 
the only satisfactory means of studying combined refraction/diffraction 
problems is to refract a wave to a breakwater zone, diffract the wave 
using an average wave speed over an arbitrary number of wave lengths, 
then refract the wave again to the shoreline. Products of refraction 
coefficients and diffraction coefficients are used to determine final 
wave heights. 

Some method of calculating combined refraction/diffraction prob- 


lems would remove these arbitrary evaluations required at present. 
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VI. CONCLUSIONS 


The basic purpose of this study was to develop a computer system 
for use by engineers in the analysis of wave behavior in coastal areas. 

The refraction program developed was based on 3-dimensional con- 
tinuous parabolic interpolation of a depth array. This procedure gives 
continuous curves through any series of data points. Adapted to 3- 
dimensions and used for refraction analysis, it appears to offer 
excellent results. 

The diffraction program is based on a form of the basic Penney- 
Price water wave diffraction methods. The system developed calculates 
diffraction coefficients and directions of wave travel which compare 
closely with manual/graphical plots. 

The entire program package offers a means of analyzing a coastal 
area in whatever depth the engineer desires in but a small fraction of 
the time required by older methods. Hopefully, this will facilitate 
more thorough design than has been economically feasible in the past. 

More testing of the programs is required; hopefully against situa- 
tions where natural conditions are known and where output could be 
compared to actual phenomena. 

Refraction/diffraction analysis is ideally suited for graphical 
plots of the output. With the present program system as a basis, output 
should be adapted to some sort of plotter or oscilloscope where solu- 
tions in progress could be watched and evaluated by an engineer. With 
console input, an engineer watching plots of his solutions would be able 
to try various design schemes for a facility and thoroughly evaluate 


their effects in but a few minutes. 
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Although improvements and further testing are needed, the system 
developed appears to offer considerable promose for better analysis of 
coastal refraction/diffraction problems than has been possible in the 


past. 
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VIII APPENDICES 


A. Fresnel Integral Calculations 


B. Program Structure 


Each program includes: a. Listing of Variables 


uP 


MAIN - 


RS RFTN 


RS INTP 


RSBETA 


RSCAFK 


RSDIFF 


RSSHOR 


b. Summary of Operations 
c. Macro-Flow Chart 
d. Program Listing 
Program to read in all commands and data and take or 
call for appropriate action. 
- Subroutine to trace a wave ray through shoaling water. 
- Subroutine that calculates depth and bottom 
derivatives at any point. 
- Subroutine which calculates ray separation factor, 
coefficient of refraction and shoaling factor. 
- Subroutine which calculates wave velocity and 
curvature at any point. 
- Subroutine which calculates coefficient of diffrac- 
tion and wave direction of travel in a diffraction 
zone. 


- Subroutine which locates and traces shoreline. 


C. Format for Command/Data Input 


D. Refraction Test on Theoretical Beach 


E. Sample Run 
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APPENDIX A 


EVALUATION OF FRESNEL INTEGRALS 


A Fresnel integral is defined by: 
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approximate C and S. For large values of t, 
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APPENDIX B 


PROGRAM STRUCTURE 


MAIN 


- Read all data and commands 

- Set constants and variables 
- Call for necessary action by 
subroutines 





RSSHOR 
| Locate and trace 
! a shorelines 


| RSCAFK 

Calculate wave speed and Ray 
curvature at any point 

(X, Y) 






RSDIFF 
Establish breakwater coordinates for use in 
calculations. Calculate coefficient of 
diffraction and direction of wave travel at 
breakwater tip or in gap. 
















RSRFTIN 
Trace path of wave ray in 
water of varying depth 









RSINTP 
Calculate depth and bottom 
derivatives at any point 

(X, Y) 


RSBETA 
Calculate ray separation factor, 
coefficient of refraction, and 

shoaling factor. 
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1. MAIN PROGRAM 


Main program to read data cards and take or call for appropriate 


action. 


PROGRAM NAME: 


Referred to as MAIN in text. 


Variables stored in Common 


AD. e e e 6 e € 


Blew, B2 . . 


BXAS, BYAS, 


BXBS, BYBS. . 


BXCS, BYCS, 


BXDS, BYDS. . . 


Ol. 6 ot 


PER «6 wt 


1) Angle in radians from +X direction to wave ray. Angle 
is measured in positive (counterclockwise) direction. 
2) In diffraction, the angle of approach of the incident 


wave ray at the breakwater. 
A in degrees; used for input/output. 


Beta, wave ray separation factor, at point preceeding 
(X,Y), at (X,Y) and at point succeeding (X,Y), 


respectively. 


Coordinates of ends of first breakwater, (A~B) on DMAT 


grid. 


Coordinates of second breakwater (C-D). If GAP is 


specified it must lie between B and C ends of breakwater. 
Wave velocity in ft/sec. 
G*T/(2*P1); deep water wave speed. 


C *T/(4*PI); referred to as CK" or K'"C in text. 


=30= 





DCDX, DCDY. 


DDDC . ., 


DMAT 6 6 e 


DXY e e e 6 


Dix, DIY . 


D2X,D2Y , D2XY 


GRID. . 


IXN,IYN . 


NBW... @ 


PL. e e 9 6 


1, 


PRRs ee le 


dC/dX, dC/d3Y in text; partial derivatives of wave speed 


to X and Y, respectively. 
dd/dc in text; derivative of depth to wave speed. 


Matrix of water depths in feet. Positive indicates water 


depth; negative is above shoreline. 
Depth at. (X,Y). (dd instext)., 
dd/dX, and dd/dY in text; slope of bottom at (X,Y). 


-d/ax’, 9°a/ae-, d-ayaxey do) cone pareiieeeons 


derivatives of bottom depth at (X,Y). 


Wave ray curvature at (X,Y). Equals inverse of radius 


of curvature. 
Acceleration due to gravity; 32.2 ft/sec. 
Size of each DMAT grid unit in feet. 


Maximum matrix dimensions being used in text. These can 
be less than or equal to dimension statement and are 
used in lieu of changing DIMENSION statement. Can use 
any DIMENSION statement which will fit in memory with 


rest of programs. 


Index for counting breakwaters; NBW=1, if none; 2 if 
one breakwater; 3 if two breakwaters (no gap); 4 if two 


breakwaters and gap between them. 
3. 1415927 
T/(4*P1); referred to as K' in text. 


(2*P1)/(G*T); referred to as K" in text. 
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Se. . © ©» © « « « Incremental distance, in grid units, between successive 


calculation points along a wave ray. 
T. «. «© © «© © » e Wave period in seconds. 


WLD .... . . j4Length of diffraction zone beyond ends of breakwaters; 
either in feet or number of wave lengths, depending on 


input. 


X,Y... + . + © Coordinates of position within DMAT which are being 


studied. 


Variables used by MAIN 


DMMY ... . - . Dummy variable, no significance. 


DX,DY,DM.. . . Used in establishing analytical depth matrices; DM is 
the depth at DMAT (IX1, IY1); DX is the increase in 
depth per X grid unit; DY is the increase in depth per 


Yeerid unit, 


fee. . -.- =. "Do Loop" indexes. 

T™X,1Y ... . . Grid coordinates on DMAT (1X,IY) array. 

ear X2 , 

Ty1,ly2 .. . .. Limits of operations being performed on DMAT. 

IWD (--) ... . Listing of decimal equivalents of "A" format characters 
for IBM 360. 

IWD1,IWD2. . .. "A" format characters read off command portion of input. 


ieee TYPl. . . . Convenient totals IX +1; TY +1. 
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LL,LL2 ... .. Index used to check through IWD(LL) array. 


LLS .... »- «. Value of LL read on last command. Needed if present 
command is a blank which indicates repeat of last 


command. 


NN... .. . . %Index used to indicate program status when calling 
RSCAFK: NN=1, calc. depth and only C; NN=2, depth is 


known, calc. © -andeeK. 


NWL. .- . - ..- © Index to indicate means of measuring length of 
ditftraction zone. (sce: Web )e 
NWL=1 if measured in wave lengths. 


NWL=2 if measured in feet. 


N1,N2,F1,F4 .. Fixed and floating point data values read off command 
card. These will be set to actual variables after command 


is decoded. 


summary of Program 


The purpose of MAIN is to read all command and data cards, then, 
either within itself or by calling subroutines, take whatever action is 
required. Its principal functions are bookkeeping, establishing variables 
for use by subroutines, and calling subroutines to take whatever action 
has been called for by a command. In some instances, variables are pre- 
determined; in some, a value is assumed unless changed by the operator; 
and in others, operator input is essential. The proper format of all 


commands is shown in Appendix C. 
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Command/Data Input 


See Appendix C for proper format. Underscore indicates Alphameric 
characters which are used in decoding command and which must occur in 
columns 1 and 10. Abbreviations which do not change underscored charac- 


ters are allowed. 


Required input for Refraction, Caleularrons. 


Either an analytical or topographic array or a combination of the two 


must be established. 


ANALYTICAL MATRIX 1X1 1X2 DX DM 


Py Ty 2 Dy 


Two cards are required. A plane surface is established which extends 
from IXl to IX2 with an increase in depth of DX per X grid unit; and ex- 
tends from LYl to IY2 with increases of DY per Y grid unit. DM is the 


depth of (1X1, IY1l). All depths are in feet. 


ieettTs OF TOPO CAE IX2 


TY1 12 
Followed by: 
TOPO IN FEET FOLLOWS 
(or) 
TOPO IN FATH FOLLOWS 
Followed by: 


Data in 10 F6.4 Format 
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This is a series of commands used in reading natural beaches into 
array which extend from [Xl to IX2 and IY1 to IY2. If topography is in 
fathoms, the program immediately converts it to feet. 


Data is read one X row at a time, 10 items per card. 


WATER LEVEL Fl 
At any time in the program, the water level may be changed to 
simulate tidal fluctuations. This card increases all water depths by Fl. 


Lower water levels can be specified by negative Fl. 


RAY DATA 1 NTC T AD xX ¥ 

This card establishes the starting point of a wave ray at (X,Y) with 
an angle AD degrees from the + X axis. NTC is any integer control number 
desired by the operator for ray identification. T is the wave period in 


seconds. 


RAY DATA 2 Bl BN 

This is an optional card which is used only if a ray separation 
factor, other than 1,00, is desired at the start of a wave ray. Bl is 
the Beta value at (X,Y), BN is the Beta value at an imaginary point 


(XN,YN) preceeding (X,Y) by a distance S along the ray. 


INCREMENT S&S 

This is an optional card which can change the initial increment 
used in tracing the wave ray from 4.0 to 2.0, 1.0 or 0.5 grid units. 
If not specified, a value of 4.0 is used. Any interval specified is a 
maximum and the program will automatically change to lesser values as 


the ray curvature increases. 
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TRACE RAY_ 

This traces wave ray specified by RAY DATA 1 and RAY DATA 2 until it 
goes off edge of array, runs onto beach, or strikes breakwater zone. 
Program will output successive X, Y coordinates of ray; wave speed, C; 
ray curvature, FK; ray separation factor, BETA; coefficient of refrac- 
tion, FKD; The Shoaling Coefficient, D; Ray Angle, AD; and the current 


increment used in tracing ray, S. 


Commands reguired if breakwaters are desired on array: 


Lack of input assumes no breakwaters; one or two breakwaters may 


be specified. 


BREAKWATER ~~ BXAS BYAS BXBS Biss 


Establishes a breakwater (A-B) from (BXAS, BYAS) to (BXBS, BYBS). 


BREAKWATER BXCS Brees BXDS BYDS 

Establishes a second breakwater (C-D) from (BXCS, BYCS) to 
(BXDS, BYDS). The first breakwater specified is automatically designated 
A-B; second is C-D. 
GAP 

Specifies a breakwater gap from (BXBS, BYBS) to (BXCS, BYCS). If 
not specified, the two breakwaters are considered independent. This card 


can only be used if two breakwaters have previously been input to program. 


SIZE GRID =IXN LYN GRID 
Sets limits of array with X from 1 to IXN and Y from 1 to IYN. If 
not specified, 50 by 40 is assumed. GRID is the size in feet of each 


grid unit of array. 
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Optional Breakwater Commands: 


WAVE LENGTHS DIFFR WLD 

This specifies the effective distance, in wave lengths, of a break- 
water beyond its tip into the diffraction zone. The wave ray will be 
traced until it either strikes a breakwater or passes within this number 
of wave lengths of a breakwater tip as measured along an orthogonal to the 
ray. Upon striking a breakwater or breakwater zone, ray values at the 


point of intersection are calculated. 


LENGTH DIFFRACTED WLD 
Similar to WAVE LENGTHS DIFFR except that WLD, in this case, is 
expressed in feet. If neither diffraction distance is specified, a value 


of 5.0 wave lengths is assumed. 


ELIMINATE BWS 
This command removes all breakwaters and gaps which have been specified. 


New breakwaters may afterwards be specified if desired. 


Input required for diffraction calculations: 


The location of one breakwater (A-B), or two breakwaters (A-B) and 


(C-D) must be specified as called for under refraction commands. 


DIFFRACT 


| > 


DIFFRACT B 


DIFFRACT C 


oO 


DIFFRACT 


DIFFRACT GAP 


Any one of these commands specifies that diffraction calculations 


should take place about the designated breakwater tip. One breakwater 
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must have already been input if A or B is specified; two if C or D is 
specified. 

If GAP is specified, an imaginary gap normal to the wave ray is 
established between points B and C. B and C must already be in the pro- 
gram. GAP also requires that a wave ray angle, A, be in the program. 

If one remains from. refraction calculations, it will be used; if not, or 
if a change is desired, a RAY DATA 1 card must be input prior to DIFFRACT 
GAP. 


Any changes in A must be followed by a new DIFFRACT GAP Command. 


DIFFR. COORD X Y 

Calculates the coefficient of diffraction and direction of wave 
travel at (X, Y). In addition to having already specified which break- 
water tip is to be diffracted through a DIFFRACT command; a wave speed, 
C; wave ray angle, A; period, T; and grid size, GRID; must already have 
been specified, GRID must be input through a SIZE GRID Command. If 
present, T and A may be left unchanged from prior operations or may be 
specified by a RAY DATA 1 card. If present, C may be left unchanged from 
a prior portion of the program or must be specified by a CALC WAVE SPEED 


or by a WAVE SPEED Command. 


Other Commands: 


WAVE SPEED C 


Sets C equal to a specified value. 


CALC WAVE SPEED XY 


Calculates and sets C equal to the wave speed at (X, Y). 


oe 








LOCATE SHORELINE 
Based on depth data, which must already have been specified, this 
command locates and traces the shoreline across an array. Shoreline 


coordinates are output. 
m(blank) & Nl N2 Fl F2 F3 F4 
Indicates repeat of preceeding command with new data values. 


ALL DONE _ 


Calls for end of the program. 
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MACRO FLOW CHART OF MAIN 


| Bookkeeping and 
isetting of constants | 


Read data card 


Check Col. 1 against allowable 
input characters (A Format) 
Go to location specified by Col. 1 


When necessary, check Col. 2 against 
allowable input characters 








Go to location specified by Col. 2 


o 





Set variables from data values on 
this or succeeding cards depending 
on command. 

Carry out command 








a 


MAIN PROGRAM TO READ IN DATA AND CALL FOR APPROPRIATE ACTION’ 


PROGRAMMED BY R D SMART AT M I T CAMBRIDGEs MASS» 1967 
DIMENSION DMAT(50+40) ,IWD(13) 
COMMON DMAT  sIXNsIYNoPI]sGsGRIDesXsYsDXYsD1LXsD1LY »sD2X sD2YsD2XxVY 9 
] AsADsTsFKsBN9Bl 2982 9CoPK s PPK »CPPK»CNsDDDC »sDCDXsDCDY » 
2 Sas NBWsBXAS sBYASsBXBS SBYPSsBXCS sBYCSSBXDSsBYDS sWLD 
199 FORMAT (A198X9A1910X%921595F1065) 
moOS8 FORMAT (10F6.4) 
197 FORMAT (19H TILT INPUT ERROR »4HWD1=9A196H WD2=39A1) 
196 FORMAT (20H DEPTH ANJUSTMENT ISsF5e3s 3HFTe) 
195 FORMAT (10H ALL DONE. ) 
194 FORMAT (18H WAVE SPEED AT X= eF5e2s 3H Y= sF5e2s 4H IS 9 
1 RO62 5M 16H FI/Ss- Car Trees sF6e394H Fle ) 
moe FORMAT (22H WAVE SPEED IS SET AT oP Gee s ohare 
192 FORMAT (12H LEST*:CAS © ssl ole PERTOD = sF 50 bs4hiseeers 
1 8H ANGLE=s F5Sels SHDEGses 2AXK=s FSecee 2H) —s oer 
191 FORMAT (20X3s2159F10-¢5) 
190 FORMAT (17H BREAKWATER FROM eF6e292H9 sF602s4H TO 1 FP Gers 
1 2Hs sF602 ) 
189 FORMAT (5H GAP ) 
188 FORMAT (24H ELIMINATED BREAKWATERS ) 
more FORMAT (26H INITIAL RAY INCREMENT IS sr Soa) 
meomrORMAT (15H DIFFR DIST IS $F Oi sOH te ) 
meme ORMAT (15H DIFFR DIST IS sF6e2s14H WAVE EENGTHS ) 
184 FORMAT (25H ANALYTICAL DATA FROM X= s1494H TO 914s 
i 8H AND Y= 914 34h 70 » 14 ) 
183 FORMAT (9H GRID IS »1T494H BY 91496H WITH SFrO67 s 
1 15H FOOT SGBARES. ) 
182 FORMAT (19H TOPO DATA FROM X= »s1l4s4H TO 914s8H AND Y= ’ 
] l4s4H TO 914) 
BPR oCFLLANEOUS BOOKKEEPING 
meet l) = —-1052753856 
TWD(2) = ~1035976640 
TWO(3) = -1019199424 
IWD(4) = ~-1002422208 
meme 5) = -—952090560 
It 6) = POPS S2545 
IWD(7) = ~-650100672 
meets) = -9185361 28 
IWD(9)= —~482328512 
IWD(10)= -431996864 
IWD(11) = —-750763968 
P9112) = -499105728 
IWD(13) = -985644992 
PT = 3461415927 
GS = 3242 
ASSUMED CONSTANTSeeceeeUNIIL RESET BY OPERATOR 
Hoe5O0 [xX = 1,50 
DO 50 IY = 1540 
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pO MATCIXslY) = =] 0000, 
NAW = ] 
BY 
B2 
a 
NWL 
WLD 
IXN 
LYN 


— 
— 
— 

—_— 


® bad Led 
® 


1 

5 
50 
40 


ent uw ow & 


READ IN ONE COMMAND/DATA CARD 


rey CY 


100 READ (559199) IWD1]1sIWD29N1°N29F15F25F35F4 
MS = LL 


CHECK IWD1 AGAINST ALLOWABLE INPUT CHARACTERSee(A FORMAT) 


CASOV TN 


DO 104 LL = 1913 

eect lWht(tt) - IWD1) 104+106;]04 
104 CONTINUE 
108 WRITE (75197) 

GO TO 138 


IF NO COMMAND IS GIVEN (BLANK)s ASSUME REPEAT OF OLD COMMAND 


i a ie a a 


Poet = LLS 


meerO LOCATION SPECIFIED BY COLe 1 


Cw OY 


MmeemGo TO (138212021375911909126:15093128513451149116514031125 142) 58 


READ A DIFFRACTION COMMAND 


Vee) 


rere o lili tlL2 = lo? 
MemclIWDO2—-IWD(LL2)) llllelll2sl1l1ll 
Prt) GO@NTINUE 
WRITE (75197) 
GO TO 138 
Sm ont | RSOIFF (LL2sIWD22Fl2F2) 
Ge, TO 100 
G 
c READ IN GRID SIZE 
C 
V2 VXN N1 
LYN N2 
GRID = Fl 
WRITE (79183) IXNs ITYNsGRID 
GO TO 100 
114 IF (IWD2 — 1077952576) 1141913691141 


c 
C READ IN GRID DATA 
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c 
1141 READ (59198) ((DMAT(CIsJ)sT=IX19IX2) sJ=lY191Y2) 
c WHEN NECESSARY CONVERT FROM FATHOMS TO FEET 
1FC(lwO2 + 1052753856) 106.0) 
fo DO 2151 1=I1X1s1X2 
Peels) J=1Y1s!1Y2 
mrol DMAT (I5J) = 6e*DMAT(I9J) 
SC uO 91 00 
Mmomel Fr (WO 2 + 1002422208) 1l6lsl47isiie) 
Meo IFC IWDO2 + 482328512) 11623116431162 


e 
C SET DIFFRACTICN IN WAVF LFNGTHS 
a 
feos WLD = Fi 
wl = | 
WRITE (72185) Fil 
GO TO 100 
‘ 
ec ADJUST WATER LEVEL 
xs 


Meo WRITE (753196) Fl 
DO 118 T=1sIXN 
Peeli8 J=15I1YN 

Pecos ATCIsJ) = DMAT(CIsJ) + Fl 
GO TO 100 


READ IN BREAKWATER COORDINATES (2 MAXe) 


irra 


MeeNGO 70 (12251242108) »sNBW 
C FIRST BREAKWATER 
122 BRAS = Fl 


BynAo= F2 
BXBS= F3 
BYBS=F4 
NBW = ? 
WRITE (75190) FlsF29F39F4 
Go, 1® 100 
SECOND BREAKWATER 
ieeoxXCS = Fil 
Beco = F?2 
BeDS = F3 
BYDS = F4 
NBW = 3 
WRITE (753190) FlsF29F39F4 
GO TO 100 


SPECIFY GAP RETWFFN 2 BRFAKWATERS 
126 TF(NBW -3) 10851273108 


127 NRW = & 
WRITE (79189) 
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GOmitGs1 00 
C 
C RFAN 
C 


IN RAY DATA 


IF(CIWD2 + 247447488) 
NTC = Nl 
= F] 
ND = F? 
AD*PI/1806¢ 


L268 
132 


r 
A 
A 
X 
nV, 
GO 


130 Bl = 
B2 = 


C 
G CALC 
ie 


ANALYTICAL MATRIX 


meet (I1WD2 + 750763968) 
moc2 1X1 = Nl 

IX2 = N2 

DX = Fl 

DM = F2 

PEeADe( 54191) 
DMAT (IX1sIY1) = 
el = {X] + ] 
TYP] = 1Y1 + 1 
Wee 1385 IX = IXP1,1X2 
DMAT (IXsIY1) = DMATCIX 
Peer387 1Y = I1YP1,1Y2 
DMAT (IX1sIY) = 
BO 13387 IX = IXP1lsI1X2 
DMAT (IXsIY) = DMAT (IX 
WRITE (79184) 
eo 710 100 


bY. Lslye2eeD Y 
DM 


v3 85 


mos 7 


C 
c ALL DONEece e QUIT 
C 
Merotmermi TE (79195) 
CALL EXIT 


CALCULATF WAVE SPFFD 


mere y C) 


137 NN = 
X 
7% 
CALL RSCAFK 
WRITE (75194) 
SO 10 100 


] 
Gl 
ie 


(NN) 


C PRESET WAVE SPEED 


NDMATCIX1s1Y - 


X9¥sCsDXY 


13091323120 


138 1s) 38203) 


— sli te 


1a 


= les Lye age 


EX ls 1 42 isles 


(FOR DIFFRACTION CALCULATIONS) 


-Q i. 








C 
fei Cc = Fi 
WRITE (791932) € 
GO 1.05100 


C 
C TRACE RAY 
C 


136 WRITE (753192) NTCsTsADoXoeY 
CALL RSRFITN (NWL) 
c Meoo NORMAL INITIAL CONDITIONS 
Bl = le 
B2 = le 
n= 4 
GO TO 100 


CHANGE NORMAL TRACE INTERVAL FOR RAY CALCULATIONS. 


wey ry 


134 S$ = Fl 
Bpoee (75187) Fi 
GeerO 100 
140 IF (TWD2-1077952576) 1401591402,1401 
meOY IF (IWD2 + 968867776) 140391 404 99403 
e 
C pee DIFFRACTION DISTANCE IN FEET 


re 
1404 NWL = 2 
mL = Fl 
MRITE (73186) Fl 
SO TO 100 
C 
fc TRACE SHORELINE. 
c 
1493 CALL RSSHOR(DMMY ) 
GO TO 100 
C 
C READ IN LIMITS OF TOPO DATA 
. 
1402 IX1 = Nl 
IX2 = N2 
ee (55191) IYIsT¥2 
Mier te (79182) IXdoeTX2s1lYisl¥2 
60 TO 190 
GC 
C ELIMINATE BREAKWATERS 
C 


142 NBW = 1 
mielhe € 73188) 
Bo TO 100 
END 
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2. SUBROUTINE RSRFTN 


Subroutine which traces a wave ray through shoaling water. 


SUBROUTINE NAME: 


RSRFIN 


Variables Used in Program: 


COMMON Variables are listed under MAIN. 


BBW1l . 


BBW2 


BBW3 


Be = . 


Pe 


DELS . 


DELTA . 


PMY . 


Dye. le 


Pp 


Ray Angle from (XN, YN) to (X, Y); an average of the 


angles at the two points. 


Last 2 preceeding values of A. Used to interpolate 


intermediate value between AN and A. 
Y intercept of Breakwater No. A-B. 

Y intercept of Breakwater No. C-D. 

Y intercept of breakwater gap. 

Y intercept of breakwater component. 
Y intercept of the ray. 

Shoaling Coefficient; D in text 


In case of intersection, the distance from (XN, YL) as 


compared to the original S$ distance. 
Change in Ray Angle, A, from (XN, YN) to (X, Y). 
Dummy; no significance. 


2*(Deep water wave length); used to estimate deepest 
reasonable depth at which waves will feel bottom. At 
depths greater than DN, it is assumed the wave ray 


curvature iS zero. 
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DN 2 e e e e 6 e 


FKBAR . . . « « 


FKN, FKNMI 
MOEN . ... « 
FLM.. . 

FLMR e eo e 
PR . . we 
FMBWl1... 


FMBW2 e e e e e 


FMBW3 e e 6 e 4 


NDXY e e e 6 e 


O.5 (deep water wave length); depth at which waves are 
substantially affected by bottom. Between DN and DN2, 
only one wave ray curvature calculation is made. At 

depths less than DN2, up to 10 iterations are used to 


arrive at a final curvature. 


Value of FK used to estimate curvature from (XN, YN) to 


(X, Y); usually it is an "average'' value. 
Coefficient of Refraction; KA Ti pexee 


The last value of FK; used to check for convergence at a 


given (X, Y) point. 


Last 2 preceeding values of FK. Used to interpolate 


average value between FKN and FK. 
Length of diffraction zone in grid units. 


Slope of the breakwater component or zone being checked 


for intersection with the ray. 

Slope of the ray from (XN, YN) to (X, Y). 
Inverse slope of FLMR. 

Slope of Breakwater No. A-B. 

Slope of Breakwater No. C-D. 

Slope (x/y) of breakwater gap. 


Number of intervals moved along ray. Ray stops at 


N=100 as a reasonable maximum. 


Index used to determine location in program prior to 
breakwater intersection test; 1 if in refraction Zone, 


2 if in deep water. 
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NFKIT «4 « « 


NP e e 6 e e 


NT e e e @ e 


SOLD Cd e e e 


pest, YDIST 


ae x2 , 


Mey? . 


Counter for the number of FK iterations at a given point. 
If not satisfactory convergence occurs within 10 itera- 


tions, an average value is used and the Fay continues. 
Intersection test index. 


Index for RSCAFK to indicate the program status when 
called; NN=l indicates need to have DXY calculated by 
RSINTP and only C calculated by RSCAFK; NN=2 indicates 


DXY is known, C and FK are desired. 


Index for RSBETA: NP=1, at N=l1; NP=2, at N greater 


than 1; NP=3, at a breakwater intersection. 


Index to assure that, when necessary, at least 2 FK . 


iterations are taken. 


Indication from MAIN as to whether WLD is in feet 


(NWL=2) or WLD is in wave lengths (NWL=1). 
Old value of S. 
X and Y components of FLEN. 


X and Y intersection coordinates of the ray with the 
breakwater component. 


Preceeding (X, Y) point. 


X and Y limits of breakwater component. 
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Functions used in Program 


BNEWE CAN. BX, CX» DX) Given 3 values, AX, BX, CX, function uses 
parabolic interpolation to calculate value of 
a point which is "DX" fraction of way from BX 


towards CX. 


coummary of Program 


RSRFIN is called by MAIN after a DMAT array is established and 
after initial ray values have been assigned, including the starting 
point, (X, Y), period, T, bearing angle, A, and the initial wave ray 
separation factors, BN and Bl. RSRFIN uses iteration procedures to cal- 
culate the path of a ray through a DMAT array of varying depths. RSINTP 
is called to calculate water depth and bottom derivatives through the 
use of 3-dimensional continuous parabolic interpolation. RSCAFK is 
called to calculate wave speed, C, and curvature, FK, at given points; 
RSBETA is called to calculate the wave ray separation factor, Beta, the 
coefficient of refraction, Rap and the shoaling coefficient, D. The 
workings of these subroutines are described in succeeding sections. 

An iteration process is used to locate (X, Y) from a preceeding 
(XN, YN) point. Point (XN, YN), the ray curvature at that point, FKN, 
and the ray angle, AN, are used to estimate a location, (X, Y), of the 
ray some incremental distance, S, ahead. RSINTP is then called to 
determine bottom depth, DXY, and the bottom derivatives at the new (X, Y). 
These data are then used by RSCAFK to calculate C and FK at (X, Y). This 
new FK is then averaged with FKN to give FKBAR. This new average curva- 
ture is then used to make a revised estimate of (X, Y). At this new 


(X, Y), the curvature is again calculated, a new average curvature is 
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determined and the (X, Y) estimate is again revised. This estimating is 
continued until FK does not change more than 0.0005 from the FK at the 
preceeding (X, Y) estimate. Once (X, Y) is established, a check is made 
to determine if the ray from (XN, YN) to (X, Y) intersected any break- 
waters or passed within any breakwater diffraction zones. If no inter- 
section occurred, RSBETA is called to calculate Beta, Ki and D. Upon 
return, the ray values at (X, Y) are printed. 

The ray is then incremented to point (X, Y) which becomes (XN, YN) 
and the process begins again. 

If at any point, DXY is negative, it is assumed a shoreline has been 
crossed. If the ray passes within 2 grid units of the array edges, the 
ray is stopped. 

In water deeper than 2 wave lengths, FK is assumed to be zero and 
Beta remains constant while the ray is incremented along a straight line. 

In depths between 0.5 and 2 wave lengths, only one calculation of 
FK is made and no iterations are used, 

As the ray is traced shoreward, S is reduced from a maximum value of 
4,0 to 2.0 as soon as FK is greater than 0.0001; again from 2.0 to 1.0 
when FK is greater than 0.0005, and finally from 1.0 to 0.5 when FK is 
greater than 0.01. S never exceeds its old value however. Initial S 
values shorter than 4.0 may be specified; this will not change the 
reduction process. 

Slope intercepts are used to determine intersection of the ray with 
any breakwater or breakwater diffraction zone (see MAIN). Breakwater 
and gap slopes and Y-intercepts are determined at the start of RSRFIN 
for use later as needed. Because diffraction beyond breakwater tips is 
considered to be orthogonal to the ray angle, these slopes and Y-inter- 


cepts must be computed individually. If intersection of a ray and 
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breakwater component occurs, ray values at the point of intersection are 
calculated directly or by using a quadratic approximation and the ray is 
stopped. 

A few safeguards are included in the program to prevent excessive 
running time in case of error. 

N is limited to 100 as a reasonable maximum number of iterations for 
a ray to require before terminating at a shoreline or an edge. 

NFK1 is limited to 10 as the maximum number of iterations allowed 
for convergence of FK values at any one point. Upon reaching 10 iterations, 
a “curvature approximation" message is printed, an average FK value is 


used, and the ray continues. 
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MACRO FLOW CHAT OF SUBROUTINE RSRFTN 


Enter from MAIN \ DXY = Deper 


L = Wave Length 
L = Deep Water Wave 
Bookkeeping, © Length 
set Constants FK = Curvature of Ray 
| S = Ray Increment 
| | C = Wave Speed 
Call RSINTP to | BW > Prealater 


Cale SID, IDi see D = Shoaling Coeff. 
K = Coeff. of Refr. 


Call RSCAFK | N 
Calc. C and 
FK 
















Cale. Beta, ee 
and other values 






DXY>0.5*L 
oe Ray 










aa Y 
: eg aries BW " = Jat intersect; out 
Sule put results 


7 Output | 
~\Results 


Close 
to edge? 





Call RSBETA 
| Calc. Beta, 


See D 


0.0005 of last 
FK estimate 


N 


N | Decrease S with 
| increase in FK 


Cale. ave ek 
(FKBAR) from 
CANIN) Comex 


Increment to next | 
Cy Ore (XN , YN) 
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peey ty) C) 


CY OVC) Camry 0) 


SUBROUTINE RSRFTN (NWL) 
SUBROUTINE WHICH WILL TRACE A WAVE RAY THROUGH SHOALING WATER. 
BASED ON 3-DIMFNSIONAL PARABOLIC INTERPOLATION. 
DIMENSION DPMAT(50340) 
COMMON DMATsIXNsIYNsPIsGsGRIDsXsYsDXYsD1XsD1lYsD2XsD2YsD2XY 
1 AsADsTsFKsBN9B1sB2sCsPKsPPKsCPPKsCNsDDDC»sDCDXsDCDY » 
2 SsNBWsBXASsBYASsBXBSsBYRSsRXCSSBYCSs BXDNS sBYDSsWLD 


BOOKKEEPING AND SETTING OF CONSTANTS. 


mencTION TO CALC INTERMEDIATE VALUE IN A SERIES 
BNFEWFCAXsBXsCXsDX) = BX+(CX-AX)¥DX/ 2041 (CX-2e*BX+AX) ¥DX¥*¥2)/26 
WRITE (739992) 
SeerORMAT (38X5 7H COEFF s- DAHSHOAa 
WRITE (75993) 
PeeoeORMAT (2X%33H N-  357H x 97H ny, > 7H A 96H e ’ 
Meo BETA »7H RFTN 35H COEF;8H CURVY GSH 0S -ee tela 
1 PK = T/(4e*PT) 
mei = 2e*P1/(G*T ) 
@e= G*T/(2-*P I) 
C = CN 
DN = 2eO*CN#T 
DN2 = O04¢5*CN*¥T 


NDXY = ] 
BN = Bl 
AN = A 
AM] = AN 
XN = X 
WN = Y 
N = ] 
NP = 1 
NT = 1 
SET UP CONSTANTS FOR BREAKWATER INTERSECTION EQUATIONS 
WATCH OUT FOR INFINITE SLOPES 
NBW = NOe OF BREAKWATERSeece 1 = NONEce ee 2 = 1 BWeeee 
2 = 2 BWeeee 4 = 2 BW + GAP 
GO TO (10510331023101)sNBW 
1901 IF (ABS(BXCS-BXBS) - 0e200001) 10015:100131002 
momma cS = BXCS + 0.0001] 
meOeeeeens = (BYCS-BYBS)/(BXCS-BXBS) 
BBW3 = BYBS - FMRW3*BXBS 
102 IF (ABS(RXNS-BXCS) - 0200091) 10035,10033;1004 
1003 BxXDS = BXDS + 0.09001 
1004 FMRW2 = (BYDS-BYCS)/(BXDS-BXCS) 
BBW2 = BYCS - FMBW2*BXCS 
103 ITF (ARBS(BXRS —- BXAS) - 0200001) 19053100551006 
mea eexAS = BXAS + 0.0001 
1006 FMBW1 = (BYBS-BYAS)/(BXBS-BXAS) 


BRW1 = BYAS - FMBW1*BXAS 
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set. ~«f 
ome 
see 





Prey Oy NAN 


a a as a sl Cre 


Us Sn 


a 


CVOV OA OY CY Ga 


FIND DEPTH AND ROTTOM SLOPES AT XsY 


mOSGALE RSINTP (DMMY) 


Sdieeree OR soHORFLINE 
RED Xe] O01] 00.) 2 
SeeckK TO SEE IF DEEP WATEResse. DEEP =e 


OS TF (DON-DXY) 12512250 


IF DEEP WATEReeNO TRIAL AND ERROR GO ON WITH FK=0sDEL BETA = O 


12 ¢C = CN 
een = FK 
FK = O« 
BN = Bl 
Bl = B2 
FKD = le 
Oe le 


CHECK FOR BRFAKWATER 


NDXY = 2 
GO TO 90 


Mom aDXY = 1 


Mees! !ONS HAVF STOPPED AT THIS POINT 
Meee OUT RESULTS 


20 AD = A * 180.0/ PI 


WRITE (79998) NoXsYsADsCsBlsFKD9DsFK 9 SsDXY 


eee ORMAT (1Xsl3e3FledsF6e2sF6e35F le 3sF oe 2s Se4sF 902s F602) 


Gimack FOR EDGFS 


ees (X—2e)*¥(FLOATCIXN) -— le = X)) 1005310053232 
Meee ((Y~2e)*(FLOATCIYN) ~ le - Y)) 1005100324 


SAFETY STOP FOR RAY AT N = 100 


eae (N=-100) 30+3051]00 


S CHECK CHANGF S VALUE DEPENDING ON CURVATURE OR RAYe 

meme Lele O2000] USE S = 4e90 GRID UNIT TSeeeTHEN USE S = 2e0siiae 
memos Toe 92005 THEN USE S = 1.460 UNTIL FK Ge Ts 96015 

ingen USE S = 9.5 

PemnOT EXCEEN OLD S VALUE 

UPON CHANGING S ADJUST B2 TO INTERMEDIATE VALUE. 


Som soOLD = S 

IF (8-32-9969) 20292023210 
202 IF (S-12-999) 20452043220 
Boer (5-0m999) 39053605230 
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210 DE tABSIFK) ~0,000]) 210.5 200.e 22) 
ie 


meOeeerTr BCABS(FK) —-. 0,005) 0222522070 
22? IF (SOLD - 22901) 30053003224 
224 S = ?e 
azo B2 = BNEWF(BN»B15A250.5) 
SOu1O. 200 
fee iF CAGBS(CFK) = Q@.01) 23232323 24C 
Mee iF {{SOLD - 146001) = O2.0001) 300.300.2234 
234 SS Vac 
be  ({ SOLD -20001) = OVO0GC 1275522257. 
mee, B82 = BNEWF (BN»B19B 260.25) 
GO TO 300 
Pa S = 0.5 
Peo —le900)]) - 0.0001) 225572255274 
Seer 6 (SOMD - 2.901) —- 060001) 23632363248 
248 B2 = BNEWF(BN>B19R2:0e125) 
c 
. TNCREMENT TO NEXT POINT BASED ON BEST FKBAR EST IMATEe 
e WHEN INITIALLY INCREMENTINGs SAVE LAST POINT AND ASSOCIATED 
cj VALUES AS A JUMPING OFF POINT FOR NEXT POINT. 
C 
300 AM] = AN 
AN = A 
XN = X 
YN = Y 
FKNM] = FKN 
FKN = FK 
NT = 1 
N = N+] 
FKBAR = FK 
C 
S Semen TO NEXT XsY POINT BASED ON BEST FKBAR ESTIMATEs 
C 
SeeepDEL TA = S * FKBAR 
A = AN + DELTA 
ABAR = AN + DELTA/?.c 
ae AN + S*COS(ABAR) 
fo= YN + S*¥SIN(ABAR) 
GO TO 10 
C 
C CALCULATE WAVE SPEED AND CURVATURE AT POINT (XY) 
a 


Bee (NTI—1) 52952980 
C Peeee FIRST ESTIMATE OF FK 
52 NFKI = 1 


wie 2 
C CALC WAVE SPEED AND CURVATURE (C AND FK) 
94 NN = ? 


CALL RSCAFK (NN) 
C N = 1 POINT ONLY REQUIRES AN FK VALUEe 
TF (N=1) 97953927045 541 
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B41 GO TO (70372582) sNT 
WO EK iowa. F K 
FKBAR = BNEWF (FKNM] »9FKNoFK 2025) 


C 
C MAKE JUST ONE PASS IN ))INTERMEDIATE)) DEPTHS 
C INTERMEDIATF = LeTe 2WlLoe OR Gel euro 
G 
IF (DXY —- DN?) 40390,90 
C 
C IN SHALLOW WATFR REVISE FK ESTIMATES JINTIL TWO AGREE WITHIN 
C O-.0005¢eeMINIMUM OF TWO FSTIMATES ARE NEEDED. 
C 
BO NT = 3 
NFKI = NFKI +1 
So 10 54 


C CHECK FOR CONVERGENCFeeeeDONT TRY MORE THAN 10 TIMES 
fee (CC ABS(FKL-FK))-0.0005) 90390584 
Serr (NFKI —- 10) 70570286 
omer ITE (73995) 
995 FORMAT (24H CURVATURE APPROXIMATED ) 
oe — {FKL + FK)/2- 


Meeck FOR INTERSFCTION WITH BREAKWATER OR BREAKWATER ZONE VOF 
OF INFLUENCEeoeeelIF THERE ARE ANY BREAKWATERSece 


CVO CON 


em (NBW - 2) 92592559255 
peor «=(ABS(X-XN) -0200001]) 2551225515256 
Meeer (ABS(Y-YN)—-0.00001 1922923256 
re CHECK FOR BREAKWATER NOsA~B 
PeorrmmMrR = (Y-YN)/(X-XN) 
BR = YN — FLMR*¥XN 
FLM = FMBW] 


ee = BBW) 
X] = BXAS 
me - BXBS 
NH = 1 

ee. 10 292 


e @eeek FOR DIFFRACTION ZONE OFF A ENDe 
257 IF (GRID-0.0001) 257453257432575 
Be ToOeWRITE (79299) 
299 FORMAT(52H PLEASE TELL ME MY GRID SIZEsOR I WONT WORK FOR YOUe) 
SO TO 100 
memroGD 10 (257192572) »NWL 
meer eeN = WLD * C * T / GRID 
@eetTO 2573 
mercer eeN = WLD / GRID 
Peete) (ABS(FLMR)-0.0001) 25783257852576 
Merem OTST = Oe 
BO TO 2577 
Ceo CMR? = ~-1¢/F WMR 
hier ot = FLEN/( 1 «+FLWR2**2 eekOe5 
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ot “Kem BXAS —~ XDIST 
YNIST = FLFN/ (1 e+FLMR#¥*#2)**® 025 
IF (BYBRS-RYAS) 26232625260 


260 YDIST = -YNIST 
mee 12 = BYAS + YDIST 
X1 = BXAS 
Y1 = BYAS 
NH = 2 
Sot O 290 


e CHECK TO SEE HOW MANY BREAKWATERS THERE AREec 
264 GO TO (92526632663272) sNBW 
¢ CHECK FOR DIFFRACTION ZONE OFF END B 


pee X2 = BXBS + XDIST 
foe bYGS — YDIST 
X1 = BXBS 
fe- BYBS 
NH= 3 
SQ TO 290 


c CHECK FOR DIFFRACTION ZONE OFF FND Ce 
268 IF (NRW-7)925979270 
Peeeeee = BXCS —- XDIST 
Bees OS -BYCS)*(BYBS-BYAS)) 270Is2702.27 02 


meoer DIST = ~YDIST 
Pee 2 = BYCS + YDIST 
mr = RPXCS 
mr = BYCS 
NH = 4 
GO TO 2909 


C CHECK FOR INTERSFCTION WITH GAP. 
272 FLM = FMBW3 


BL = BBW3 
Xl = BXBS 
X2 = BXCS 
NH = 5 

mee tO 292 


C Memeo kK FOR INTERSECTION WITH BREAKWATER NOeC=D 
274 FLM = FMBW2 


ee = BBW2 
X1 = BXCS 
X2 = BXDS 
NH = 6 

Bo TO 792 


C meek FOR DIFFRACTION ZONE OF END De 


eee = BXDS + XDIST 
fm = BYDS —- YOIST 
X1 = BXDS 
mime - BYDS 
NH = 7 
C Meee SLOPE INSERSECTION (M*#X + B) TO SEE IF RAY AND 
c BREAKWATER COMPONENT INTERSECT. 


meer (ABS(X2-X1) -0.00001) 2901 3299232902 
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2901 XP? = XP? + O2AN0] 


BIO? Filip = (Y2-Y)])/(X2=%)) 
Bl = 1 > FLM*X] 
Boel r {ABS(FLMR=-FLM) - 0.00001) 2921572922522 
eo7) MeMR = FlMR + 0.0001 
2922 XI = (BL-BR)/(FLMR-FLM) 


tet A AM XI—XN) ) 29452935235 
Qe ((X2—-X1)*(XI—X1)) 294592965296 
294 GO TO (25732649268 927452745216592) sNH 


C 
G Meee Y DO INTERSECT CALCULATE POINT OF INTERSECTIONS Aaa 
. INCLUDING RAY ANGLE» RAY SEPARATION FACTORs 
. COEFFICIFNT OF REFRACTION AND SHOALING FACTORe 
= 
meee = (XI—XN)¥CY-YN)/(X-XN) + YN 
DELS = ((XI-XN)¥*¥2 + (YI-YN)¥*2)*#0.5/5 


B1 = BNEWF(BN»sR19B2eDFLS) 

A = BNEWF(AM1 5ANSASNELS) 

coe A @ 180./PI 

ED (1-/B1)* #025 

X I 

ve I 

NN = ] 

CALL RSCAFKI(NN) 

IF (NP=1) 296252952:2961 
Z96G2,D = le 

SO 10 2963 
7961 NP = 3 

@atlt RSBETA (NPsD>FKD) 
meow RITE (75298) XIisYI 

Mier eE (73297) 

WRITE (75993) 

WRITE (75998) NoXeVYsADeCeBloFKDsDoFK 9S sDXY 
298 FORMAT (41H INTERSECTION WITH BREAKWATER ZONE AT X= ’ 

1 F6e3s4H Y= si cie > ) 
297 FORMAT (36H RAY VALUES AT INTERSECTION FOLLOW. ) 

Se 1O 100 

92 GO TO (971519) sNDXY 


a 
=< << il 


CALCULATF BETA 


ry erenee 


peoeeKN = FK 
Peet RSBETA (NP»D>sFKD) 
Boe tO 20 


MO OY 


END OF RAY CALCULATIONS- 
POOUWRITE (75994) 
SOGFORMAT (13H RAY STOPPED ) 
RFTURN 
END 
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3. SUBROUTINE RSINTP 


Subroutine which calculates depth and partial derivatives of the bottom. 


SUBROUTINE NAME: 


RS INTP 


Variables used in Program 


COMMON variables are listed under MAIN 


DMMY . 6 9 


Be, JY x 


Dex. DLY 


fee TY. .« 


NX, NY. . 


IXNY, IYNY 


many) . 


DD (NX) . 


DD1 (NX). 


DD2 (NX). 


é 


Dummy variables, no significance 

Integer X and Y; not rounded. 

Decimal fraction of (X-JX) and (Y-JY). 
(JX-2) and (JY-2). 

Counters for X and Y. 

X and Y grid points. 

Bottom depth at location DMAT(IXNY, IYNY). 
Bottom depth at (IXNX, Y). 

od/oY at (1 xipen 


a*d/dy* at (1XNX, Y). 


Functions used in Program 


DEPTHF(Al, A2, A3, A4, DELTA) 


Calculates data value at point of interpolation. Given. 
four equally spaced data values, Al, A2, A3, and A4, 

this function uses continuous parabolic interpolation, 

as described in the text, to calculate a data value at 
point DELTA. DELTA is that fraction of the distance from 


point A2 towards A3 at which we desire data values. 
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DEP1F (Al, A2, A3, A4, DELTA) 
Calculates slope of curve at point of interpolation 


using continuous parabolic interpolation. 


DEP2F (Al, A2, A3, A4, DELTA) 
Calculates second derivative of curve at point of 


interpolation using continuous parabolic interpolation. 


Summary of Program 


Continuous parabolic interpolation is reduced to three functions; one 
to calculate a data value; one to calculate a slope; and one to calculate 
a second derivative. All values are at any abcissa value between the 
center two of four equally space data points. 

With (X,Y) located within the center grid of a 4 x 4 data array, 
these functions are used by 4 interpolating parabolas, parallel to the Y 
axis, to define 4 points on a plane parallel to the X axis and passing 
through (X, Y). Using these 4 new points, and their derivatives when 

parabola 
necessary, the same functions are then used by an interpolating/parallel 
to the X axis, and passing through (X, Y) to completely define (X, Y) 


and all desired derivatives of a surface through (X, Y). For a detailed 


explanation of this process see Section II-E (C) of text. 
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MACRO FLOW CHART OF SUBROUTINE RSINTP 


f Enter from RSRFTN Es 
\ or RSCAFK 


{| Misc. Bookkeeping | 
| and setting constants} 


Calculate values and derivatives 
| for 4 "Y intercepts" 





Calculate values and derivatives 
| for (X, Y) using 4 "Y intercept" 
values 






=i 





merry © Y CY CY) 


VAOMAN 


CVNet) 


CVO) 


SUBROUTINE RSINTP (DMMY) 


SUBROUTINE CALLED BY RSRFTN OR RSCAFK WHICH WILL USE 
3-DIMENSIONAL PARABOLIC INTERPOLATION TO CALCULATE 
DEPTHsBOTTOM SLOPESsAND OTHER PARTIAL DERIVATIVES 
OF ROTTOM. 


DIMENSION DMAT (50940) .9(4)+DD(4) 50114) = Dpee 

COMMON DMAT sIXNsITYNoPI2GsaGRID 9X 9 Y¥ sDXY sD1X sDLY 9D2X s9D2YsD2XY 9 
1 AsADs TsFX »BN9Bl 9B29CaPK sPPK »CPPK »CNsDDDC»sDCDX sDCDY » 

2 S»NBWePXAS »sBYASS BXBS eo BYRS®9BXCS »BYCS 9 BXDS sBYDS sWLD 


BOOKKEEPING AND SETTING OF CONSTANTS. 


FUNCTIONS FOLLOW 
DEPTHF (AL 2A22sA39A49DELTA) = A2+(DELTA/2¢)*(A3—-Al)—-(DELTA#*2/2.-) 
i ® (A446 ¥A3B +50 FA2—-2e%*AL)+(DELTA¥#3/2 6) ® (AG—3 0 *¥A34+30%A2-Al ) 
DEP1F(A12A29A32A49DELTAI= (A3-Al1)/2e-(NELTA/2 ©) ¥(A4—-5 0 FAB +70 *A2 
1 —3e*A1l) + (DELTA**%2) * (A4—-3¢*¥A3+3e*A2-Al ) 
DEP2F(A1 9A2 9A3 sA4sDELTA) =A3—-2 e*#A2Z+A1+DELTA* (A4—-3 0 *A3+30*A2-Al ) 
JY=Y 
BEY = Y-(FLOAT(JY) } 


JX = X 

Dex = -Xo>- (FLOAT(IX)) 
IX = JX -2 

Ries JY — 2 


Mee RPOLATE 4 TIMES PARALLEL TO Y AXIS FIRSTececses 
MeECULATE DEPTHS AND DERIVATIVES FOR Y-INTERCEPTS OF Shee 


DO 200 NX=194 
IXNX = IX+NX 
Poms O NY=1354 
mrNY = JY + NY 


180 D(NY) = DMATCIXNXSIYNY) 


DOWUNX) = DEPTHF(D(1)> Dl2)> Dl3a)>s DO(4)> DEY) 
DN1(NX)= DEPIF (D(1)29 Dl2)% Dl(3)9 N(4)s DLY) 


Peep D2°NX)= DEP2F (Dll1)s Dl2)s DI(3)>5 Dl4)> DLY) 


Meee ESULTS OF LAST 4 INTERPS TO CALC VALUES AT (Asiv® 
DXY = DEPTHF(DD(1) »DD(2) sDD(3) »DD(4) sDLX) 
OLY = OEPTHF (DN1(1)99D1(2)2DN1(03)sDN1(4) »DLX) 
Pe “SRP IPP (OD2Z(1)sDD2(2) »>DD2(3)sND2(4) 5DEX} 
Pix=  DEPITF (DD(1) »PD(2)sDD(3) »DD(4) »DLXx)} 
fe. =sperze (DN (1) »~DD(2)+D0(3).00(4) sDLA) 
D2XY=DEP1F (DD1(1)sDD102) »0D1(3)»DD1(4) »DLX) 
400 RETURN 


END 
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4, SUBROUTINE RSBETA 


Subroutine which calculates the ray separation factor (BETA), the 


coefficient of refraction, Rae and the shoaling coefficient, D, for a 


ray. 


SUBROUTINE NAME: RSBETA 


Variables used in Program: 


COMMON Variables are lised under MAIN. 


D2CDX2 


BZCDY 2 


D2CDXY 


D2DDC2 


FKD . 


PED . 


PP, QQ 


Wave group velocity. 

Deep water wave group velocity. 
Shoaling coefficient; D in text. 

3° C/3X- in text. 

97C/ay~ tn text. 

97C/3X9Y in text, 

aaa dee in text. 

Coefficient of Refraction; Ka itiee sts 


Convenient summation. 


p and q of equation for B2 shown under summary of program. 


Summary of Program 


After the ray trace subroutine, RSRFTIN, has determined an (X, Y) 


coordinate, RSBETA is called to use a "Finite Difference" method, (See 
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text, Section I-E (d)), for projecting from the two most recent Beta 
values to the next succeeding value. 

The program requires 2 initial Beta values to start, one at the 
present (X, Y) point, Bl, and one at the preceeding (XN, YN) point, BN. 
All points must be equally spaced. 

Each time called, RSBETA increments itself from the prior BN, Bl 
values and calculates the succeeding B2 using 

= ese re psa 
oa 2+p 5) a 5 + p s on 

It must be kept in mind that this method is projecting an 
estimated Beta value one step ahead of ray tracing. One return from 
RSBETA, Bl is the present (X, Y) Beta value. 

In addition, the program calculates the Shoaling Coefficient at a 


point by using 


_con, 1/7 


CG? 





where 


ao a/ Lb 


sinh (4 T a/L)? 


CG =0.5* C* (1 + 


and CGN is the deep water value of CG. 


It also uses the Beta value to calculate the coefficient of 


refraction at a point. 
= 
/ 


K, = (Beta) 2 


d 


maeprogram notation this is: 


FKD = .g1) 72/4 


See; — 





MACRO-FLOW CHART OF SUBROUTINE RSBETA 


(Enter from RSRFIN 


Yes 





1 No 
| Bookkeeping and increment 


from last Beta values 


if 












| Calculate necessary partial 
i derivatives and use them to 
find peand q. (Pe andsoo) 











Solve = a8 ps-2 
B2= em 2+ps 2 a+ SZ) 


find EKD, (K ,) 
| FKD=(B1) + 


1/2 





CG=0.5¢C1: sinh (4rd/L) 


Yes 





D=(CGN/CG)1/2 


CGN-CG 


=e — 





SUBROUTINE RSBETA (NP D>sFKD) 

CG SUBROUTINE TO CALCULATE RAY SEPARATION FACTOR (BETA) FOR 
DIMENSION DMAT(50,40) 
COMMON DMATsIXNoITYNsPIsGsGRID 9X s Ys DXY sD1LX eDLY sD2X sD2Y xs D2XY 5 
1 AsADsTsFK»BNBl 9B29CoPK »PPK »CPPK sCNsDDDC»sDCDXsDCDYs 
2 S»yNBWesBXAS sBYASs BXBSSBYRS »sBXCS sBYCS 9 BXOSSBYDS eWLD 


INCREMENT BETA TO PRESENT POINT 


erty () 


GO TO (700,700,705) »NP 
700. BN Bl 
Bl ez 


ott 


STILL IN RELATIVELY DEEP WATEReeeBETA = 1¢000000 


wre 


IF (ABS(CPPK—-1.20)-0.209000001) 70457023702 
704 B2 = le 
GO TO 710 


CALCULATE NECESSARY DERIVATIVES 


AAD 


PK*(2e*(PPK/(1.2+CPPK) #RIPK/ (1 seGPPK)}) + Gi PPR aa 
Pe=—CPPKT**2 — PPR®*¥2/) «+GRP ie A aazee) 

D2X/DDDC - (D200C2*D1xX**2 ) /DDDGHF3 

D2Y/DDOC = (D20DC 2D | Yen 207 DOLG- 

DZAY/7DDDOC — (D2DDC2*0) Xap] ¥) 7 PPG =5 


eeize D2DDC2 
] 
D2CDX2 
BZEDy 2 
PZGDXY 


tou ua tt 


SOLVF BETA EQUATION 


CORAM) 


fate Gree OS( A) + DCDY*SIN(A) 7G 

QQ (D2CDX2*(SIN(CA)*#*¥2)— (D2CDXY*#2e¢*SIN( A) *COS(A)) 

] tee eC Yoh (COS AI) wre?) Ae 

BZ = ((PP*S-22)/(PP*¥S+2-))*BN + ((4e=20*QO0OKS*42? ) /( PPS Tze | 
SRE RETURN PRINT OUT Bl AS BETA FOR POINTS Asy¥ 


CALULATE COEFFICIFNT OF REFRACTION. CPKRe 


OOO) 


FKD = (14¢/R1)**0.5 


MmeeGULATE SHOALING COEFFICIENT. (0) 


Vea 


705 FPD = 4e*PI*DXY/(C#T) 
mom —emmear ec #(1+(FPD/(0.5*( EXPCFPD)—EXP(—FPD) J) 
GO TO (70637085708) sNP 

706 CGN = CG 


D = le 

NP = 2 

SO 10 710 
708 D = (CGN/CG)**0.5 
Pero RETURN 

END 
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5. SUBROUTINE RSCAFK 


Subroutine which calculates wave velocity and ray curvature at any 


given (X, Y) point. 


SUBROUTINE NAME: RSCAFK 


Variables Used in Program 


COMMON variables are listed under MAIN. 


CL ..... . Wave speed used to check convergence of iterations. 
LL ..... . Index; no significance. 
NN ... . . « Index to indicate program status. If NN=1, call 


RSINTP to determine DXY, then only calculate C. If 


NN=2, use present DXY and calculate both C and FK. 


Summary of Program 


When called with NN=1, RSCAFK calls RSINTP to determine water depth, 


then determines wave speed using an iteration process to solve: 


peo 271d 
C On tanh L 


The program then returns. 
When called with NN=2, depth data have already been determined. 


Using these data, it determines wave speed then goes on to determine ray 


curvature using 


ac ac 
(sin A AG - cos A (sy)? 


ale 


FK = 


=) — 





The program then returns. 

As a safeguard, a maximum of 80 iterations are allowed for C to 
reach successive values within 0.0005 of each other; otherwise a "Wave 
Speed Approximation" message is printed and the program continues using 


an average C value. 
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MACRO FLOW CHAT OF SUBROUTINE: RSCAFK 


iEnter from MAIN | 
or RSRFTITN | 





Call RSINTP to 
—| calc. depth at 



















Calc. ray curvature 
it oC aC 
FK= C (sin AG) cos A (y)? 
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Cree res formrercye yt y OY CY 


Pye y 


Caren 


VO OY 


SUBROUTINE RSCAFK (NN) 


SUBROUTINE WHICH CALCULATES WAVE VELOCITY AND RAY CURVATURE 
AT ANY GIVEN (XsY) POINT 
USUALLY CALLED BYSRSRFIN DURING RAY TikAGeeeeerPUIATIONS. Sear 
MAY BE CAINBER BY MAIN JUST TO CAGE WAVE VEROC iy 
AT ANY SPBGr TED LOCATION. 


DIMENSION DMAT (50:40) 

COMMON DMATsIIZANe PYNsPI2GesGRID»sXsYsDXY silxXs Diy sDZ24A3D271sD2Zx"> 
1 AsADsTsFKsBN2Bl »B29C2PKsPPKsCPPKsCNsDDDC»DCDXsDCDY » 

2 S»oNBWeBXAS sBYASSBXBS SBYBS SRBXCSSBYCS 3S BXDSsBYDSsWLD 


merCALLED BY rNeeneeD TO HAVE DEPTH CALCULATED FiksSI. 


GO TO (53355) »NN 

53 CALL RSINTP (DMMY) 
GN = G@ *°T Vette Pl!) 
Cc = CN 


CALCULATE WAVE SPEED9sCe 


comet = C 
DO 56 LL = ITsee 
© = CN * TANM@Q2er-ey XY) / CrCl )) 
IP CABS(CL=C) 20.0005) 58558356 
pemct = (C + CLI72Z. 
WRITE (723999) 


999 FORMAT (26H WAVE SPFFD APPROXIMATEDe ) 


IF CALLED BY MAINeeeGO BACK TO MAIN 


Seman TO (59557) eNN 


SPeEGULATE CURVATURE OF RAY (= FK ) 
Peer PK = C ™% PPK 
Pome = (PR) 8O( 262 CPPK)7(1e-CPPK2)+ALOG (1etCPPK) 
1 = AoOGa (ls marr Ns ).) 
DCDX = DIX/DDDC 
DCDY = DIY/DNDC 
mx = (Cle/C)I*¥(CSIN(CA))*DCDOX — (COS(A))*DCDY) 
Bees | URN 
END 





6. SUBROUTINE RSDIFF 


Subroutine to calculate Coefficient of Diffraction and Direction of 


Wave Front Travel at end of semi-infinite breakwaters or at a breakwater 


Gap. 


SUBROUTINE NAME; 


RSDIFF 


Variables used by Program: 


COMMON variables are listed under MAIN, 


ee, BB. . 
AAl, BBl . 
ABW e 6 © ® 
ABWL ¢ é ¢ 
AP ® ® e ¢ 
Al, A2 

B1X, B2X . 
BXA, thru 
BYD ® ¢ 6 ¢ 


A, B in text. 
Used to save "half" solution across breakwater gap. 


Bearing of breakwater or gap, in radians, relative to 


X axis, 
ABW value used for direction checks. 


Bearing from breakwater tip or center of gap to (xX, Y) 


and relative to. asic. 


Interim summations. 


For diffraction at one end of breakwater, the desired 

end is assigned coordinates, (BXB, BYB); the other end 

is assigned coordinates (BXA, BYA). For diffraction at 

a gap, imaginary coordinates BXA through BYD are assigned 
based on an imaginary gap normal to the wave ray (see 


text) 
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EC, saWicea ea. 
i 
CIM, CIP —; 

DADTH, DADR , 
DADX, DADZ . 


DADZ1, DADXL, 


DBDZ1, DBDX1 


DAl, DA2, 


Dei, DB2 .. 
DBDTH, DBDD . 
DBDX, DBDZ ., 


DFM, DFP, 


moe etc. . . 


DFMDR, DFPDR 


DFMDTH , 


PeePTH . . « 
MEE 5 5 8s 
1 
ED gl 


MeGZR ., ., 


C in text, part of Fresnel Integral. 
Coefficient of diffraction. 

cos (TH-THN), cos (TH+THN) respectively. 
dA/d08 and dA/dr in text, 


dA/OX, OA/OZ in text. 


Used to save "half" solution across breakwater gap. 


Interim summations. 
0B/d8 and OB/oOr in text. 


oB/dX, dB/oZ in text. 


Appropriate derivatives used to make passes through 


similar equations without rewriting equations. 


OFM/or, OFP/dor in text. 


OFM/08, OFP/08 in text. 

Gap length. 

Angle of wave travel in radians. 
Angle of wave travel in degrees. 


dg'/dr; also do/dr; see Note. 
Note: Same calculations are needed with both o and o',. 
A first pass is made through calculations (M=1) using 


g values although the variable name indicates o' values. 


ee 





DSIG2T 


Br. 


DU2DR, 


DW2DR 


DU2DTH 


DW2DTH 


Dx. DY 


Dae DYT 


FACFP, 


FBCFB 


FACFPX 


FBCFPX 


FACFPZ 


FBCFPZ 


FF1, FF2 


FLTK 


FM, FP 


6 


Upon completion of calculations o, results are saved, oa' 
values are inserted, and a second pass is made (M=2). 

First pass, with o, calculates Ul, Wl and their deriva- 
tives; second pass, with o', calculates U2, W2 and their 


derivatives. 
d56'/38; also 00/00; see Note following DSIG2R. 


Half width of imaginary gap normal to wave ray. 


dU2/dor, dOW2/dr or dUl/dr, dOW1l/dr. See note under 


DSIG2R. 


0U2/098, dW2/08 or 3U1/08, W1/d0. 


Various X and Y differences used locally in various parts 


of program; no significance. 


Local subtotal; no significanre. 


cos (FP), sin (FP) in text. 


d(cos(FP))/dX, d(sin(FP))/dX in text. 


d(cos(FP))/dZ, d(sin(FP))/dZ in text. 
Interim summations, 
K in text; 2 *PI/(C*T). 


FM, FP in text; kr cos (TH-THN), kr cos (TN+THM) 


respectively. 
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FN. 


FT. 


el, 2 


IWD2 


Ca 


NGAP 


NGAPS 


NN, 


NQ . 


n in text, index for Fresnel summations. 


; 2 
t in text; To ly 


Variables from MAIN; coordinate locations (X, Y). 


Alphameric variable from MAIN; the end of breakwater 


which is being diffracted, A, B, C, or D. 


Index passed from MAIN which indicates desired action by 
subvoutame., If LL=1l, 2, 3, or 4$0d#itrrwetion is at Ase, 
C, or D end of breakwater respectively. If LL=5, diffrac- 
tion is at Gap. If LL=6 or 7, perform diffraction for 


given coordinates, 
Index; see Note following DSIG2R. 


Index (1, 2, or 3) used to return to correct location in 


program after calculating a bearing. 


NGAP=1 for diffraction at the end of a breakwater. 
NGAP=2 for gap diffraction with higher values being 
assigned during the process of gap solution in order to 


keep track of the status of calculations. 
Original value of NGAP to be restored at end of program. 


NN=1 if imaginary gap is being established at start of 
program; NN=2 at end of program when gap is being restored 


for future use. 


Index to indicate whether angle AP is greater than 


angle A (NQ=2), or less (NQ=1). 


r in text; distance from (X, Y) to either a breakwater 


tip or center of gap. 
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ot 


S14 


SIG] 


SIG2 


SS 


TERM 


Te 


THN 


THD, THND 


THE L 


ie U2 , 


Wl, W2 


one YC 


ZZ 


ne 


coor o'. 


G Ine rexe, 
o' dn@texe. 
S in text. Part of Fresnel Integral. 


Terms of summation series. 


Angle from breakwater to (X, Y) as measured at breakwater 
tip. 
Angle from breakwater to ray of incident wave as 


measured at tip. 
TH and THN in degrees. 
Or (8+6,) or 02> (8-8) 3 no significance. 


Used to save "half" solution across breakwater gap. 


Same as text. 
Coordinates of center of gap. 
0.5 (2*R/(C*T)); no significance. 


Convenient summation. Used in sin and cos approximations 


of Fresnel Integrals. 


Zs — 





summary of Program 


Prior to starting diffraction calculations, incident wave data must 
be input to the program as previously described in the section on MAIN 
commands, 

In order to solve a diffraction problem, RSDIFF is first called by 
MAIN after a DIFFRACT Command has been read. Depending on whether told to 
diffract breakwater and A, B, C, D, or Gap, index LL is set to 1-5 
respectively by MAIN. Then RSDIFF is called to assign breakwater co- 
ordinates (BXB, BYB) to the end of a breakwater at which diffraction is to 
take place and assigns BXA, BYA to the other end. BXAS through BYDS are 
never disturbed and are kept for later use if desired. The main portion 
of RSDIFF is such that it always computes diffraction about the B end of 
an A~B breakwater. When called to diffract a gap, the program first 
determines an imaginary gap, (see Text Secttton III, B) which is normal 
to the incident wave ray and sets index NGAP=2 to indicate a gap problem. 
Control then returns to the main program. 

When later called by MAIN after a DIFFR COORD command has been read, 
RSDIFF calculates the coefficient of diffraction and direction of wave 
travel at (X, Y) and about the previously designated tip or gap. 

The main portion of the RSDIFF program is arranged to calculate 
values needed to solve a standard semi-infinite breakwater problem. 
Breakwater gap solutions can be built out of semi-infinite solutions. To 
solve for a gap, it is therefore necessary to make multiple passes through 
the body of the program with "semi-infinite breakwater" portions of a gap 
then add the results for the final desired values. 

Equations for this program are explained in detail in the text of 


this report. Due to their length and complexity, they are not repeated 
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in the summary which follows; but to the maximum extent possible, similar 
notation is used in the programs and in the text to facilitate cross 
reference, 

Using text notation TH, THN, the program first determines necessary 
angles, and R. These are then used to find FM, FP and their theta and R 
derivatives. 

Next SIGMA 1 and SIGMA 2 are derived and used in Fresnel integrals to 
determine C and S. These are then used to give Ul, Wl, U2, W2 and the 
derivatives of each with respect to both Theta and R. 

At this point, one has all the input needed to find A, B, and their 
Theta and R derivatives. 

Continuing the building process, the derivatives of A and B are used 
to find 0A/dX, dA/dZ, JB/dOX, and dB/9dZ. 

If the solution is only desired for a semi-infinite breakwater, 
these values are used directly to find the coefficient of diffraction and 
the wave direction of travel. 

If a gap solution is required, A, B, and their derivatives are saved 
as a partial solution and coordinates (BXB, BYB) are established to coin- 
cide with the "center" of the breakwater gap. It is then necessary to go 
back and repeat those steps necessary to solve for FP and its derivatives. 
This information is used to determine exponential derivatives needed for 
the solution. 

Coordinates (BXB, BYB), (BXA, BYA) are then changed to coincide 
with the C-D breakwater. All steps are repeated as required to find A, 
B, and all their X, Z derivatives again. 

It is then possible to add all partial solutions for the gap to give 
totals for A, B, and their derivatives. These are used to find the 


coefficient of diffraction and the wave direction of travel. Upon 


7 = 





completion of a gap solution, coordinates of the imaginary breakwater are 


restored so that further calculations can be made if desired. 
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MACRO-FLOW CHART OF RSD 


Enter from MAIN ) 


IFF 


Go to location specified 
by second character of command | 







A,B,C,D 






Set up (BXB,BYB) 
and (BXA,BYA) to 
perform diffrac] 
tion off "B" end 
of breakwater. 







Establish imaginary gap 
normal to wave ray 

Set up diffraction off 
iB” end, 








Set up coords to| 
domcaligne art 
center of gap 











Calc. SIGMA2 
and its 
derivatives 






Restore Gap 
if used 


Output 
Results 


er > 


Save Ul, | 
W1, etc. 














Determine ] 
IAP ,ABW,R 
Determine 
TH , THN 


| Cale. FM, FP 
and their 
derivatives 


ente 
of gap _ 
? 
No 


Calc. SIGMAl and 
its derivatives 


Evaluate Fresnel Integrals 
to get C and $ 


Galc. U2, W2 and 
their derivatives 











Save values _ 
Set up to 
do calc. 

"GO ences 





Or 









Ye 
SIGMA] 
Values 
Calc. K' and Wave 51a 2 
Dir. of Travel No 
A Yes 


Add up gap | No 'B 
values r end 
9 
Ves 
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Calc. A, B, and their 
derivatives (Two passes | 


through same equations) 





Calc. DADX, DBDX, | 
DADZ, DBDZ 





ame Cy Caneel ) 


Crier 0) CY 


G 


G 


SUBROUTINE RSDBFF (bGeRWD2.F1>5F2) 


PROGRAM TO CALCULATE COEFFICIENT OF DIFFRACTION AND DIRECTION OF 
WAVE (FRONT) TRAVEL AT END OF BREAKWATERS OR AT A 
BREAKWATER GAPe 
DIMENSION DMAT(50940) 
COMMON DMATsIXNoIYNePI sGesGRIDsXsYasDXYsD1X9D1Y9D2XsD2YsD2XY5 
1 AsADsTsFKsBN5Bl1sB2 9C s9PK »9PPK sCPPK sCN yDDDC sDCDX sDCDY » 
2 SsNBWsBXASsBYAS »BXBSsBYBS sBXCSsBYCS sBXDSsBYDS sWLD 
GO TO LOCATION SPECIFIED BY COMMAND . 
GO TO (500155002 3500325004 35005 5500635006) sLL 
BOOKKEEPING TO SET UP DATA FOR CALCULATIONS 
PROGRAM CALCULATES VALUES OFF B END OF BREAKWATER 
FOR SOLVING GAP» THREE PASSES WILL BE NEEDED TO SOLVE 
FOR THREE IMAGINARY ENDS « 
DIFFRACTION AT A END OF BREAKWATER 
5001 BXB = BXAS 
BXA = BXBS 
BYB = BYAS 
BYA = BYBS 
NGAP = 1 
GO TO 5008 
C DIFFRACTION AT B END OF BREAKWATER 
5002 BXB = BXBS 
BxA = BXAS 
BYB = BYBS 
BYA = BYAS 
NGAP = 1 
GO TO 5008 
DIFFRACTION AT C END OF BREAKWATER 
5003 BXB = BXCS 
BYB = BYCS 
BXA = BXDS 
BYA = BYDS 
NGAP = 1 
GO TO 5008 
DIFFRACTION AT D END OF BREAKWATER 
5004 BXB = BXDS 
BYB = BYDS 
BXA = BXCS 
BYA = BYCS 
NGAP = 1 
5008 WRITE (735099) IWD2 
5099 FORMAT (30H DIFFRACTION CALCULATIONS FOR 4 Als 
1 19H END OF BREAKWATER.e ) 
RETURN 
DIFFRACTION IN BREA“WATER GAP 


C 


“SO 





c 


‘a 
. 


C 
G 
C 


c 
C 
€ 


€ 


ESTABLISH AN IMAGINARY GAP. 


5005 


5009 


5007 


20 15 
2 O18 


5097 
5999 


Does 


NGAP = 2 
NN 1 


p 
DX = BXCS ~- BXBS 


DY BGS" —e ov GS 

DG (DX*¥*¥2 + DY*¥*¥2)**005 
mC (BxXSS.* BXGS)/ 2 
ie (Byes) + BYCS) 72 e 

N = 3 

GO TO 203 

DT = DG/2e * SIN(A-ABW) 
DXT = DT * SINCA) 

DYT = DT * COS(A) 

BXB = XC = DXT 

BXA = XC = 3e*DXT 

BXC = XC + DXT 

BXD = XC + 3e*DXT 

mee = VC + DYT 

Seeke = WE a SeetDYT 

orc = YC = DYT 

BEyDs= YC “wee *DYT 


GO TO (501595016) 4NN 
WRITE (755098) 
FORMAT (46 H DIFFRACTION CALCULATIONS FOR BREAKWATER GAPe  ) 
WRITE (795097) 
FORMAT (52H BASED ON IMAGINARY GAP ASSUMED MORMAL TO WAVE RAYe) 
WRITE (755999) BXBsBYBsBXCsBYC 
FORMAT (21H IMAGINARY GAP FROM) 9F60e391H99F6e39 
4H TO »9F6e391Hs9F603) 
RETURN 


SET COORD POINTS 


5006 


X Fl 
io EZ 
NGAPS = NGAP 


DETERMINE BEARING AND DISTANCE FROM BWB TO XsY_~ (APo9R) 


900 


903 
5011 


N = 1 
lo Say = Bene 
wR = KR- XS 


R = (DX*¥¥2 + DY**2)**0e5 * GRID + 020000001 
We (ABSUDXOS0900001) SOMSO92 +5012 

ABW = PI/2-e 

GO TO 5014 


5012 ABW = ATAN(DY/DX) 
QUADRANT CHECK 
59014 IF (DY*DX) 50455045506 


yar 





=Set@ it a? ° 





EAC? MAN NAN NAN 


aa 


MANDAN 


SOS TECDX) 507% s50sa508 
506 IF(DY) 50795093509 
507 ABW = ABW + PI 
Ge 10 5098 
508 ABW = ABW + 2e*PI 
PF ( ABW2 oGPI) 509%502 5502 
502 ABW = ABW —- 2e*PI 
209) GO TO (51 1s5 3s 5007 15N 


DETERMINE BREAKWATER BEARING (ABW) 
511 AP = ABW 
DY = BYA = BYB 
DX = BXA = BXB 
N = 2 
GO TO 503 
DETERMINE ANGLE FROM BW TO RAY(THN) 


513 THN = ABS(A=ABW) 
DETERMINE ANGLE FROM BW TO XoY (TH) 


TH = ABS(AP=ABW) 
mT) 51695169515 
515 TH = 2e*PI = TH 
owe CPAN=P 1578 9518 s517 
af ThNe= 2etPT = THN 
SET MISC SUMMATIONS TO BE USED LATER 
318 FLTIK = 2e*PI/(C#¥T) 
ZS LETRA YR/PI YK F065 
CTM COS: TTH-THN) 
cre SOG { | HASKIN ) 


CALC FMs FP» AND THEIR DERIVATIVES e 


FM = FLTK¥R*¥CTM 
FP = FLTK*R#¥CTP 

DFMDR = FLTK*CTM 

DFPDR = FLTK*CTP 

DFMDTH = =FLTK*R*SIN (TH=THN) 
DFPDTH = -FLTK*R*¥SIN (TH+THN) 


NGAP = leeeeNO BWe99 NGAP = 2eeeeDOING LEFT SIDE 
NGAP = 3 eeweeDOING CENTER9999 NGAP = 4 wee eDOING RIGHT SIDE 
GO TO (519951995775519) sNGAP 


USE SIGMA VALUES TO CALC Us Ws AND THEIR DERIVATIVES 
THROUGH THE USE OF FRESNEL ‘INTEGRALS. 
SEC SIGIMPVALUES FIRSTweeelmanN SIG2Z™VAEUES “(USE LOOP) 
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NAN 


NAN 


519 SIG2 = 2225s Net On oz TH<THN)?> 
THX = 0«65*(TH-THN) 
he OSIG2) 5 eer, Sao 5 92 
5191 NO = 1 
Ge TO 5193 
a192 MQ = 2 
59193 M = 1 
SO PT = PISSIGZHAZ/ 26 
RE (ABSCTPAK) — 08990001) 5101951025 5802 
mt DSTG2T = 22 
GO TO (510495105) 9M 
meee DSIG2T = =-ZZ 
GO TO 5104 
peeoe DSGIG2ZT = G®GZ*COS( THX) /(2e*SIN( THX) ) 
9104 SI = ABS(SIG2) 
EVALUATE FRESNEL INTEGRALS (CC AND SS) 
TP (S1Siye) “S125540 ,540 
ioe iaele Se USE SERIES SOLUTION 
See ol = STP /2.)**¥0-5) 
S14 = SI*#4 
FN = le 
heEkM = SI 
Cer= Sl 
514 FN = FN+1le 
TERM = — TERM*¥S1I4%(4e*FN—70)/ (16 0*FNEX3—52 oe KFNEKX2454.%FN—-18. ) 
Ge= CCHTERM 
Tim CABSMORERM) - 0.0001) 5245)45514 
pee = CCX((20/PI )**0-5) 
FN = le 
SS = (SI**3)/36 
TERM = SS 
SoU FN. 2 FNYT e 
TERM = -TERM*SI4*( Ge FNS ee )/ (160 kFNXX3B—28 ge KFNXX2 +14 e%XFN—2e ) 
SS = SS+TERM 
Mm (ABS (TERM)=-0-0001) 53835309530 
538 SS = SS¥*((2e/PI)**005) 
GO TO 545 
IF SI GeTe 3¢ USE SIN/COS APPROXIMATIONS 
540 Z2 = (PI/2e)*(ST*¥2-26/(PI*SI)*%*¥2) 
GGe= 0b * SIN» CG2y/(PI*Sdo 
eo = O«> = COS Bz) 7tP lS!) 
CALC U2sW2,AND ALL THEIR DERIVATIVES. 
545 U2 = Oc S*(10—SS—CG) 
W2 = 00e5%*(SS-CC) 
IF (ABS(SIG2)-02000001) 54519545235452 


9452 DSIG2R 
DShGe2 T 


(SIG2/(2e%#R))*(SIG2/(ABS(SIG2))) 
DSIG2T * (SIG2/(ABS(SIG2))) 
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AOAOAD 


ADAODDNNA 


C 


5451 
5453 


GO TO 5453 

DSIG2R = SIG2/(2e#R) 

DU2DR = -O45*(SIN(FT)*DSIG2R + COS(FT)*DSIG2R) 
DW2DR =O04e5*(SIN (FT) *DSIG2R-COS (FT)*DSIG2R) 
DU2DTH 20% 5 HS I Ne FT EDS IG2T + COS (EI? *DS1G2T ) 
DW2DTH OeS=GSIN (FY) *DSIG2ZT - COS Saiieiaes G27) 


tt ol 


ON FIRST PASS RESET VALUES TO SHOW Ul oW1sAND THEIR DERIVATIVES 
SET SIGMA2 AND GO BACK TO CALCULATE U2 9W2sAND THEIR DERIVATIVES. 


0 


GO TO (5503560) »M 
Ul = U2 

Wl = W2 

DU1DR = DU2DR 
DW1DR DW2DR 
DU1DTH = DU2DTH 
DW1DTH = DW2DTH 
STG2e= —20eB2S WN ( OwS* ETHETHN Je) 
THX = 065*(TH+THN) 
M= 2 

GO TO 510 


USE SAME BASIC A AND B EQUATIONS TO CALCULATE BOTH R AND THETA 


at 


560 


we 
i 


I 


1 


1 


DERIVATIVES OF A AND Be 


RST DO R DEVIVATIVESe 

M=1] 

DFM = DFMDR 

DFP = DFPDR 

DU1 = DUIDR 

DU2 = DU2DR 

DW1 = DW1DR 

DW2 = DW2DR 

Fae = Ul * COSTPM) +°W1 *® SINCFM) 

A2 = U2 *COS(FP) + W2 *SIN(FP) 

Bl1X = Wl *COS(FM) = Ul * SIN(FM) 

oem = We * COStTEP) —- U2 * SINUFP) 

DA1 = DUI *COS(FM) = UL*¥SINCFM)*DFM + DW1*SIN(FM) 

+ W1*COS(FM)*DFM 
Brie = DUI2 * COSHFP) = U2*SIN(FP)*DFP + DW2 * SINMFP) 
ree COS (FP Rep rr 

Pel) — bw) * COSCEM) = W1*SIN(FM)*DEM — DU1*SIN( FM) 
=—U1*COS (FM) «DFM 

DB2 = OW2*COS(FP) = W2*SIN(CFP)*DFP - DU2*SIN(FP) 


= "U2"COS (FP YeDFP 


WATCH OUT FOR PROPER QUADRANT 


570 


ate 


GO TO (5703580) sNQ 
DADTH = DA1 + DA2 

DBDTH = BB1l + DB2 

eer TO (57125575) »M 

AA = Al +* A2 
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Bee BX 4 B2x 


SEM) O <555 
SOO IRPENGAP=1) 58 Oar O!s5611] 
5801 DADTH = =SIN(FM)*DFM = DAI] + DA2 
DBDTH = =COS(FM)*DFM - DBl1 + DB2 
G@ 710 3812 


5811 DADTH = DA2 -DAI1 
DBDTH = DB2 -DBl 
pem2 GO TO (5825575) 5M 
Dare “RR UNGAP=] ) SRSZ Ih'82 1 95822 


5821 AA = COS(FM) - Al + A2 
BB = =-SIN(FM) = B1X+ B2x 
ae TO 265 
5822 AA = -Al + A2 
BB = =-B1X + B2x 
= 
C SAVE R VALUESeeeeSET VARIABLES SO AS TO CALCULATE THETA 
c DERIVATIVES AND DO IT OVER AGAINe 
G 
So Me= 2 
DADR = DADTH 
DBDR = DBDTH 
DFM = DFMDTH 
DFP = DFPDTH 
DU1 = DUILDTH 
DU2 = DU2DTH 
DW1 = DW1DTH 
DW2 = DW2DTH 
GO 10 365 


e 
e CALC DERIVATIVES OF A AND B WITH REGARDS TO X AND Ze 
c 
975 DADX = DADR*¥SIN (TH) + DADTH*COS (TH) 7 R 
BeDX“="RGURESIN (TH) + DBDOTH*COS (TH)/R 
. WATCH OUT FOR QUADRANTS 
9s GO TO (57519575295751:95751) »NGAP 


ore CRDZ = —DADR*COSTPHY + DADTH*SIN(TH)/R 
BBDZ = —-DBDR*COS(TH) + DBDTH*SIN(TH)/R 
eer lo 54 
5751 DADZ = DADR * COS(TH) - DADTH*SIN(TH)/R 
DBDZ = DBDR*¥COS (TH)-DBDTH*SIN (TH) /R 
C 


C IF THIS IS A GAP GP BACK AND DO AGAIN FOR CENTER AND RIGHT SIDE 
C 
5754 GO TO (574957655775578) sNGAP 


¢ 
c PevE OLD "VAUUES FOR LEFT SIDE OF Bw 
C 
976 DADZ1 = DADZ 
WS071 = WepZ 
DADX1 = DADX 





DBDX1 = DBDX 


TH11 = TH 
AAl = AA 
BB1l = BB 
DO CENTER NEXT 
NGAP = 3 
BYA = BYB 
BXA = BXB 
BXB = (BXB+BXC)/2e 
BYB = (BYB*BYC)/2. 
GO 710 500 
C 
es SAVE CENTER VALUES (ONLY A PARTIAL SET) 
C 
577 FACFP = COS(FP) 
fab CF Pe =e Sanlit FP ) 
Bri #2 DRPPDORFSINCTH) + DFPDTH * COS(TH)/R 
Pe2 = DPPDReCOStTH) - DFPOTH*SIN(TH)/R 
FACFPX = =-SIN(FP)*FF1 
mocre PX = COSKEP ) FF | 
ACE Ze = 9S BFP ) FF 2 
meee PZ = =COS(FPIAFF 2 
S GO BACK AND DO RT SIDE OF BW 
BXA = BXD 
BYA = BYD 
BXB = BXC 
BYB = BYC 
NGAP = 4 
GO TO 500 


IF THIS WERE A GAP PROPLEM MAKE SUMMATIONS.- 


DOES POINT LAY BETWEEN GAP TIPS OR BEHIND A BREAKWATER. 
MAKE BW SUMMATIONS 
eee OTH A -— PT)*{(TH = Ad) 57995799581 


AN OM 


581 DADX = DADX + FACFPX 
DBDX = DBDX + FBCFPX 
DADZ = DADZ + FACFPZ 
UpeZ = DBDZ + FBCFPZ 


AA = AA + FACFP 
BB = BB + FBCFP 


379 DADX = DADX1 + DADX 
DADZ = DADZ1 + DADZ 
DBDX = DBDX1 + DBDX 
DBDZ = DBDZ1 + DBZ 


MADE TTceeeeceeWHEWecccccce 
NOW CALC DIRe OF TRAVEL AND COEFFe OF DIFFRACTION 


-VSaQ@e 


CY OY Cras 





/otrre 





< 
C 


OY OM 


OO 


574 IF (ABS(AA*DBDZ-BB*DADZ )-0200001) 574125742 95742 
5741 DIR = PI/26 
GO TO 5744 
5742 DIR = ATAN((AA*DBDX-BB*DADX) / (AA*DBDZ-BB*EADZ) ) 
DIR IS WRT BWeeeeeRETURN TO XsY COORD SYSTEM 
CORRECT FOR NEGATIVE ATANS IN SECOND QUADRANT. 
5734 IF (DIR) 574395744 55744 
5743 DIR = DIR + PI 
5744 ABW1 = ABW 
IF (ABS(A=ABW1)-PI) 57489574895749 
5749 ABW1 = ABW1-2-e*PI 
5748 IF ((2e*PI-(A-ABW1) )*(A~ABW1) )57463574695745 
5745 DIR = DIR + ABW 
GO TO 5747 
5746 DIR = ABW = DIR 
cee IF «6(DIR—2e%P1) 57319573295732 
eee DIR = DIR=—2e*P I 


Siem DIRD = DIR*180,/P' 
GDIF = (AA@™2 + BReFEZ)*EO,5 
PRINT OUT RESULTS 
BRITE (755096) 
PRITE (7935095) ADseTsC 
WRITE (795094) XoY 
WRITE (735093) DIRDsCDIF 
NGAP = NGAPS 
PeonmoORE ARTIFICIAL.GAP IF USED 
NN = 2 
IF (NGAP=2) 50162500935009 
5096 FORMAT (30H WAVE CHARACTERISTICS  ) 
5095 FORMAT (9H ANGLE 9F7e2s 18H DEGes PERIOD s9F 602s 
1 Headlpot Coa. SPEEDs FSO2> 9H Tiieeees,.) 
5C94 FORMAT (14H AT POINT X = 9 F5e2s SH Y= 9 F5e2) 
59093 FORMAT (27H DIR SOR TRAVEL Ts sF7e2s 
1 29H DEGesCOEF OF DIFFRACTION IS mee fant) 
RETURN 
END 
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7, SUBROUTINE RSSHOR 


Subroutine which calculates shoreline location from a depth array. 


“ 


SUBROUTINE NAME: RSSHOR 


Variables used in Program: 


COMMON variables are listed under MAIN. 


D..« .«. ec. . ##=*Fraction of distance between two grid points at which 


shoreline crossing has been interpolated. 


DX, DY «.. .. D distance along X or Y grid line respectively. 
INC. « «4. c© e Index used to keep track of which of 6 possible grid 


crossings is being checked; INC=1 during original search 


for a shoreline segment. 


IX, IY .. . . Coordinates at primary point of shoreline search; 
program looks for shoreline between (IX, IY) and 


Cixi, IY) cr between (IX, ye and (ix ty 499 


1X2, IY2 .. - Coordinates at secondary grid point opposite (IX, IY). 


Shoreline check occurs between two points. 


NER({LL).... A check list of coordinate locations at which grid 
crossings have already been found. Coordinates are 
stored as an integer (100 * IX + IY); with a plus value 
indicating grid crossing between (IX, IY) and (IX+ 1, IY), 
and a negative value indicating grid crossing was between 


(IX, IY) and (IX, [y+ 1). 
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NN Do loop counter, no significance. 


NNS Index; NNS=l1 if last grid crossing was along Y grid line; 


NNS=2 if last grid crossing was along X grid line. 


NS Index; NSs+1l if searching between IX and IX +1; NS= -l 


if searching between LY and IY +1. 


NSL Index of total number of grid crossing located and 


Stemededm, NCK(-—)e 
XD, YD Coordinates of shoreline at a grid crossing point. 


so. YS First coordinate points along a shoreline segment; used 
to close loop if shoreline traces back to origin as on 


an island or enclosed bay. 


Summary of Program 


When called by MAIN, RSSHOR finds and traces shoreline segments 
across a DMAT array by using the fact that, only at a shoreline, the 
product of two adjacent depths (DMAT (1X, IY)) * (DMAT (1X2, IY2)) is 
negative, 

Starting at (2,2), the product of adjacent grid depths is determined 
along grid lines parallel to the Y axis until an initial starting point 
of a shoreline segment is found. Once a shoreline is found, it is then 
traced to its ending, either at the matrix edge or back at its start such 
as at an island. 

Each time a shoreline grid crossing is found, the coordinates of 
(IX, IY) are stored as integer (100 * IX + IY) in NCK (NSL). It is given 
a negative value if the grid crossing were between (IX, IY) and 


(IX + 1, IY); it is given a positive value if the crossing were between 
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(IX, IY) and (1X, IY + 1). Each possible crossing point is first checked 
against NCK to be sure that it is in fact a new crossing point, not a 
duplicate, 

Shoreline traces are made by checking the six adjacent grid bounds 
for a new crossing; when found, it is stored in NCK. A check is then made 
of its 6 adjacent grid bounds. The NCK listing eliminates any line backing 


up upon itself. 


Each time a shoreline is ended, a new search starts at (1,2) until 


all segments have been found. 


Output is a series listing of consecutive shoreline coordinates. 
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MACRO-FLOW CHART OF RSSHOR 


Enter from MAIN 
Bookkeeping a) 


Search for increment. 
of Shoreline 


INC= 1 
TX= 2 
















No Yes 


Return 





, Bdge — pe TX=1x42 


= eae 
Tvs? 
TY2=1v4+1 Ty=Tytl 







DMAT (IX, IY) 
*DMAT (IX2, 1Y2) 





Eliminate 
zero depth 










(IX,IY) already 
on NCK list? 











Trace Shoreline 
Check adjoining 
grids bounds 

(6 possible) 





Interpolate to find 
Dist to crossing 
polite) 









Go to X or 
Y search 
(NNS=2) 

X or 
crossing: 
(NS = 1,-1) 














Print crossing 
point 






DY=D SS 

Store (IX,IY) in 
NCK (+ if along Y 

erid, = adreavene X 


grid 
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ere, (Yty OVO OM 


CY 


mw) ay) OO) arr). CVOX 


DNOODHN 


CY OYA) 


SUBROUTINE RSSHOR (DMMY ) 


SUBROUTINE -wWHkGHmWhbet LOCATE THE SH@RELINE 
BASIC APPROACHeeePRODUCT OF ADJACENT DEPTHS IS NEGATIVE 
WHEN CROSSING SHORELINE) 
USE STRAIGHT LINE INTERPOLATION BETWEEN GRID POINTS OF APPOSITE 
SIGN TO APPROXIMATE SHORELINE 
WSE NG@K eeeexX) AS CHECK LIST FOR COQORDINWATE LOCATIONS 
BOOKKEEPING 
DIMENSION DMAT (50940) 5NCK(200) 
COMMON DMATsIXNsIYNoPIsGsGRIDsXsYsDXYsD1XsD1LYsD2XsD2Y sD2XY 9 
1 AsADsTsFK9BN5B1lsB2 9CaPK sPPK »9CPPK sCNsDDDC sDCDX sDCDY » 
2 SsNBWsBXASsBYAS 2 BXBSsBYBSsBXCS sBYCS 9 BXDSs8YDSsWLD 
NSL = @ 
ZER® OUM GAECK LIST 
DO 402 LL = 1.4200 
402 NCK(LL) = O 
SEARCH FOR INCREMENT OF SHORELINE. 
INC = 1 ee NOT IN PROCESS OF TRACING SHORELINE SEGMENT 


404 INC = 1 


IX = 1] 
496 NS = 1 

IX2 = IX 

ey ie 2 


410 I1Y2 = IY+l 
BBIGe CHECK 


420 IF ((IX"1)*CIXN-IX)) 49094229422 
eer CCT YSKLPECILYN=1Y)) 49094259425 


CHECK PRODUCT 
425 IF (DMAT(CIXsTY)*DMAT(IX2s1Y2)) 430942653411 


RerereemeNnT TO NEXT "POINT 
feer BNCREMENTINGe UNT bm ENTIRE GRID IS CHECKEDeec.s 


411 IF (INC=1) 46094125460 
412 ITY = IY+l1 

IF (TYN=ITY) 41494143410 
414 IX = IX + 2 

IF (IXN=-IX) 49534069406 


SeriTiiNnnatesANY ZERO DEPTHS 
426 IF (DMAT(IXsIY)) 42894279428 


427 DMAT(IXsITY) = DMAT(IXsIY) + 060001 
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‘3. 


‘=F 





ryt. 


OY 


NAAN 


NAN AN 


CY NON 


NANDONM 


428 IF (DMAT(IX291Y2)) 42594299425 
429 DMAT(IX2s1Y2) = DMAT(IX291Y2) + 920001 
GO TO 425 


GRECK FILE FORSPREVIOUS LISTING OF Geiaieig 
430 DO 445 NN = Io»NSL 


IF (NCK(NN) = NS¥(1O0O*IX+1TY)) 44594119445 
445 CONTINUE 


INTERPOLATE BETWEEN GRID POINTS 
WAS IT AN X OR Y CROSSINGeee 
D = DMAT(CIXsIY)/S(DMATCIXsIY) - DMATCIX29s1Y2)) 
IF (NS=1) 450345259452 
450 DX = D 
DY = Oc 
NNS = 2 
GO TO@s5s 
452 DY = D 
DX = Ovc 
NNS = 1 
STORE GRID POINT IN FILE 
srreome + VALUE IF ALONG Y GRID 
STORE = VALUE IF ALONG X GRID 
455 NCK(NSL) = NS¥*¥(1LOO*¥IX + IY) 


XD = FLOAT(IX) + DX 
YD = FWOAT(TY) + Dy 
IF (INC#1) 457594565457 
456 XS = XD 
> ee YD 
WRITE (79496) 
496 FORMAT (5Xs22H SHORELINE COORDINATES) 


PiReeNT CROSSING POINT 


457 WRITE (79497) NStL»xDsYD 
497 FORMAT (15917X92F8.3) 


NSL = NSL+1 
INC = 2 WHEN FOLLOWING SHORELINE 
INC = 2 


SEARCH FOR ADJOINING SHORELINE GRID CROSSING 
L@@K IN 6 POSSIBLE DIRECTIONS 


460 GO TO (4593470) »NNS 
459 GO TO (462 3462 5463 5464 5465 946634675468) s INC 
G62 1X = [X¥] 
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INC = 3 
461 IX2 = IX 

NS = oi 

GO TO 420 


463 IX = [X-2 
IF (1X) @65254652 54631 
4631 INC = 4 


GO TO 461 
464 INC = 5 
ITY¥2 = IY 
469 IX2 = IX+] 
458 NS = -1 
GO TO 420 
4652 IY2 = IY 
465 INC = 6 
IX = IX+1 
GO TO 469 
466 INC = 7 
ye 1 YY 
Ty2 = IY 
GO TO 458 
467 INC = 8 
IX = [Xl 


IF (IX) 46859468 9469 
DID WE GO BACK TO STARTING POINTeee 


468 IF (ABS(XD-XS)-2e) 481594819490 
481 IF (ABS(YD~YS)=2¢) 48234823490 
IF BACK AT STARTs PRINT STARTING COORD. 
482 WRITE (79497) NSLsXSsYS 
GO TO 490 
470 GO TO (4725947294173 54745475 947694775468) 5 INC 
472 TY = IY+l 


we = 3 
Ome ly¥2 = IY 
NS = =-1 
GO TO 420 
famemlyY = IY-2 
INC = 4 
GO TO 471 
474 INC = 5 
IX2 = IX 
Om? 1Y2 = I¥#1 
480 NS = 1 
GO TO 420 
475 INC = 6 
LY see 
GO TO 479 
476 INC = 7 
IX = IX+] 


Sy ara = 





IX2 = IX 


GO TO 480 
477 INC = 8 
ey = LY md 
GO TO 479 
C 
C END OF SHORELINE SEGMENT 
C 


490 WRITE (79499) 
499 FORMAT (5Xs 22H END SHORELINE SEGMENT  ) 


G GO LOOK FOR ANOTHER SEGMENT 
GO TO 404 
495 RETURN 
END 
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APPENDIX (C) 
FORMAT FOR COMMAND/DATA INPUT 


CARD COLUMNS 


i= 20 Cee ake 6 Oe | — 41) 41-50 leo 0 61-70 
ANALYTICAL MATRIX Ext 1X2 DX DM 
ee. lY2 ony 
feels OF TOPO x TX2 
Texel by2 


meee IN FEET FOLLOWS 
mero IN FATH FOLLOWS 


INCREMENT = 

RAY DATA NTC ii AD X " 
RAY DATA 2 Bl BN 

WeeTER LEVEL Fa). 

menace RAY 


Oe te SHORELINE 


BREAKWATFR BXAS BYAS BXBS BBS 
BREAKWATER BXCS Bes BXDS BS 
GAP 

See ORI D IXN LYN GRID 

Pee ENGHTS DIFFR WLD 

Pom DIFFRACTED WLD 


ELIMINATF BRKWTRS 


teats RAC T 
PPE ERACT 
DIFFRACT 
Pipe RACT 
DIFFRACT 


Groat) a > 


AP 


Pier Rk COORD X i 


(CONTINUED) 


are 








CARD COLUMNS 
20) Oe 2625053 1-—40 


i a 
PAG wWAVE SPEED 














ANALYTICAL 
Meotl CASE 


OM NANA NP WNMNFeH ZZ 


X 
3eN00 
Bee ao 
7224) 
8.2654 

10.066 
1124/76 
122883 
14.2285 
is 08 | 
172966 
mi. (53 
18.2436 
192114 
mo. (35 
207118 
202776 
ene 102 
212424 
216/744 
224060 
O26! 2 
22468] 
222984 
Gere 25 3 
2.5/6 
232864 
24.144 
24.2417 
244682 
24.938 
252184 
Pen 4 19 
254640 
254847 
Zoe 0 36 
26.2205 
264349 
262458 


REFRACTION TEST ON ANALYTICAL BEACH 


DATA FROM X= 
1 PERTOD = 


" 
2,000 
ie CZ 
7244 
BeGo? 

10.076 
126 Se 
144242 
LS e272 
Lee oe 
176743 
18.674 
19.2409 
20215] 
202524 
Cnc 9 
Ces! fT 
212656 
222938 
220423 
2224810 
226 200 
2320594 
28a 99) 
244392 
242/97 
250206 
AD Oe 0 
262039 
264463 
2ow O92 
eel 2s 
Zee OF 
Cee 16 
Ziti 6 7 > 
ae 46 
292696 
Ciera elots 
302573 


fer STOPPED 


A 
45.090 
G5 6 Onk9 
452040 
45.074 
452125 
45622? 
452369 
454589 
45.2917 
46402 
462/28 
472115 
472579 
48.134 
48.2452 
482/798 
oe ks 
492586 
BOOS. 
5045290 
512050 
514627 
5202535 
D249 3 
eueoo > 
54.2495 
Dost go 
560346 
572403 
20% 002 
599-484] 
614253 
624824 
642585 
bbe 85 
68.886 
712610 
146994 
79.2709 
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(D) 


SEE | ABLE 


palo 50 AND Y= 1 TO 40 
ta Oi E Ce ANGLE= 452¢0DEGeX= 3-200Y= 
COBRMmOHOAL 
C BETA RFTN @ocr ~ CURVY > 

20049 12-000 124000 1200 02-0001 4-200 
202e¢48 1-009 1-000 1-00 02-0002 4.00 
Amo el «6 00} 1eCOO0 14400 O2CCO3 2-00 
20046 1400}1 046999 1200 02-0004 2-00 
20643 12002 Og999e0 me 0007 2.60 
20040 12-004 02998 0299 04209010 2-200 
20¢35 14006 04997 04698 040016 2-00 
20427 12-010 O8e95 0.928 (Wee 2c 
20615 126016 80m O2 0.97 BO eee ee. oD 
19697 14024 GRBs 04696 OPOUSie 
19685 12029 SOeeeG 0.96 Os00G wie 
19e71 146036 "Oe C.95 Us00 74 ae 
19254 14043 046979 04294 O0O-20089 1.200 
196.33 1.052 | ORS 30-94 0.01906 12:00 
19621 1.058 “Oepe 0.94 0.0115 O.e0 
19.08 12¢0673 (Wie 0.9: 0-0126 Us5C 
160693 1¢069- Ose. 73) «0607 Ue 
16078 1.076) Oe eee os 0.0150 9.50 
18660 1083 “OC S680. c2 _°G20163 0.50 
18e41 1209] (ORGS 7 30.92 O.0lt7y Ose0 
18.21 12099 Uae e ol “0.0 1o2et.so 
17298 JelO8  Osee O50. 92 Dse210 0259 
Llel3 Vell? Cees 92 O60229 0.c0 
17e¢45 14128 Ones 0691 0.0249 0.650 
L7elS 14139 Oey. 9) ) «€6OSOC2 TL Meee 
16682 14150 0.932 9.91 0.0296 0.50 
16045 124162 Oue2z7 9292 O.032370.-0 
16605 Jel?6&. O6e22 0692 O5045200.c0 
15660 124191 046916 0692 0.40386 0.250 
P5011 12206 O91] O694 590.6474 70. 
14655  1le222 0690500 ee 4oc os 
13293 14239 O4898 0094 040519 0.50 
T3623 10257 O88 C2 ee 0s O57To 0550 
l2etsele2 (6 One eeeRS (060552 0250 
Lleol oie 297) Ostrom 7060745 0650 
10643515318 OseeeeO4 O40869 0.50 
eres 1 341 02854 1210 02-1048 O¢50 
7e00 10365 Ombeeeie20 0.41347 0650 

De FS ] 2390 0-848 1042 Oe2045 0250 
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3 


O90 


Oe eH 
562000 
502 34 
G7e5] 
4468 
41.85 
39-01 
2 Ora lb 7k 
Baeg 2 
30a ae 
Zia De 
26611 
24265 
rae le 
21270 
20295 
20620 
19245 
18669 
Lie 
1715 
1638 
L260) 
14-8] 
14292 
Leie2 2 
1264] 
ies? 
10276 

9eoZ 

9.207 

Be22 

1034 

6046 

50 56 

4065 

Seah 3 

2019 

le 65 

(Vee 








APPENDIX 


al) 


SAMPLE PROGRAM RUN 


DIFFRACTION OPERATIONS CALCULATE VALUES FOR AN ASSORTMENT 


MeeoCELLANEOUS CASES-e 


REFRACTION OPERATIONS SHOW VIRGINIA 


INPUT DATA DECK eccocccves 


Maye SPEED 

Ze GRID 50 
BW 

Bw 

RAY DATA 1 2 
DIFFRACT A 

Siete COORD 


RAY DATA 1 1 
DIFFR COORD 


Pi FRACT PD 
DIFFR COORD 


RAY DATA 1 i 
DIF FRACT GAP 
Cite COORD 


ELIMIN BREAKWATERS 

BW 

DIFFRACT A 

RAY DATA 1 1 
Poet COORD 


Pilea ACT B 
DIFFR COORD 


FLIMINATF BRKWTRS 


20% 

BO. 
106 
2260 
Ge 


106 
12e 
“Le 
6 e 
605 
Te 


l2e 


Dis 
266 
Cle 
2705 
286 
40 


2le 
236 
2065 


106 
909 
106 
1061 
7 


20 6 
200) 


Sig— 


BEACH iz sisS 


106 
226 
G0 


146 
136 
14e 
1356 
Ize 
125 
136 
146 
136 


306 
Bie 
336 
3305 
34 6 
90. 


266 
2 Oke 
236 


106 


SO) 
los 
156 
156 


156 
15 
13% 


VoEe FiGe 


206 
306 
le 


206 


OF 


4) 


206 
306 
lie 


10s 


Ze 








few. vD ee Be 





S25 GRID 50 40 320006 
ANALYTICAL MATRIX il 50 O« = 
a3 40 O« 
mel TS OF. TOPO ] 50 
1 32 
TOMO IN FATH FOLLOWS 
95 95 96 97 97 98 98 99 96 95 
92 90 8 8 B5 8 3 82 82 82 82 82 
82 82 82 82 Bz 8] 82 82 eZ 82 
8 3 B3 84 Bi5 85 86 86 87 87 88 
88 89 98 92 92 93 93 94 93 92 
92 93 94 95 95 96 96 97 96 94 
92 89 87 84 8] 80 80 80 Bel Bell 
8] 8 ] 81 80 80 80 80 8O 8] ome 
(INTERIM CARDS OMITTED) 
—45 meee () -40 -40 =40 -40 -40 -40 a =-40 
7G) = 400C«C“—C CCRC 0 SO HAO « =40 eo 
=O “40 =a6 -40 —40 =a = 2h. =e soe = 1b 
BMS) 44 LS 45 fee? aS 458 46 465 Lot 
=o -40 a) -40 -40 =e -40 2p} -40 =—20) 
-40 —~40) —-49 -40 -40 -40 -40 —-40 =A.0 ot) 
0) —4f) oa -30 =—29 = 20 -15 8. =-00 09 
12 19 22 28 30 on 36 38 40 405 
G5 4] 45 42 425 Pees 5 44 445 “5 
ac) —-49 —-40 -40 —-40 (59 -40 Sent) =aitE) =O) 
=2.0 -40 -40 -40 -40 ee -40 ee -40 -40 
=4.0 —40 -40 -40 -40 -40 ~35 =) =25 2) 
Za) 108 -C0 08 12 19 20 Ze 25 29 
CUNTEXREMeCARDS OVITTED) 
San) ate) ero -40 -40 -40 -49 — 1256 — 20) -40 
-40 eds () =e -40 —-40 —40 -49) =a) =—£:) =: 0 
_ =a = &() =e) -49 -40 —49% -4° —-40 -40 = ce 
met) —t[) a) -4O0 —-40 —4% -40 = =40 -40 
-40 -40 ee -49 -40 —40 —40 = e6() ay) =—26 
may eOATA 1 a Le 3206 500 
TRACE RAY 
eee DATA 1 a one ae 620 
TRACE RAY 
RAY DATA 1 Be Le 3206 720 
TRACE RAY 
omer SHORELINE 
fee DONE 
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OUTPUT INFORMATION DECK eocccces 


web SPFED 15 SFT MAT 20.<005RT/7S 
Gea Da 1S 5 @eRY 40 WITH 80-90 FOOT SQUARES. 
RREAKWATFEFR FROM 1@500, 10.00 TO 200Gtieee 20200 
BREAKWATER FROM 22.005 22.209 TO 30.005 30.200 
DIFFRACTION CALCULATIONS FOR A END OF BREAKWATER. 
WAVE CHARACTERISTICS 
ANGLE 90.00 DFGes PERTOO 4s600 bees SPEERZOCOG Fil7 sae. 
AT POINT X = 8e00 Y= 14-09 
DIRSG@raTRAVEtL «IS 92037 DEGse GC OLFPFRACT TON 13 ier eo 
WAV EPECHAR ACTER TST TCS 
ANGLF 90.00 DFGes Pak TOD 4200 SFCes SPEED2Z0. 00 iy = aa. 
DR wore Tie. TS 82.60 DEGss COEF OFM Ie ecTlON To Gee]. 
VeVESChHAinaGIreRITSTICS 
ANGLE 90.00 DFGes» PeRlOD igQOr SE Ces SreeUzc0s00 Fl /SeGe 
AT POINT X =12290 Y= 14.209 
DIRSCe MRAVEL TS 6267 DEGssCOEF OF DIFFRACTION. a2 
WAVE CHARACTERISTICS 
ANGLE 135200 DEGe» PERIOD 4200 SECes SPEEDZ0.00 “Fly sae. 
AT POINT X = 6600 Ye 12.009 
DIR®OF TRAVEL IS 13/641 —DEGS.COEF OF DIFFRACTION, 12 iw: 
WAVE CHARACTERISTICS 
Bore 1354600 DEGes PERIOD the SE Ces SPEEDZOS 008 I aa. 
AT POINT X = 6650 Y= 12.59 
DOR Or IRAVE IS 1346020 DeGesCOEF- OF DIFFRACTION 1S) Once 
weve NCHARaGTERISTICS 
Prramee 135.00 DEGe>» PERIOD 4.00 SECes SPEEDZ0 00 2F 7 ane. 
AT POINT X = 7200 Y= 13-400 
BIRVOF TRAVEL 15 128280 DEGssCOEF OF DIERRACTIOND] ee. 
WAVE CHARACTERISTICS 
myete 135.00 DEGe> PPREOP 4 OOS ECs SPEED20.00 FIZ are. 
AT POINT X = 82-00 Y= 14.00 
PIR OFSIRAVEt: IS 115600 DEG.s COEF OF DIFFRAC] LONG 30.2 
WAVE CHARACTERISTICS 
Poss 125.00 DEGsw> PERIOD 4 «00S Rees SPEEOCUS COs Fr I7SEae 
AT POINT X =12200 Y= 13409 
PROF TRAVELS IS 56631 DEGe sC@EE OFS DIFPPREAGCIION 1S” Oclzea 
DIFFRACTION CALCULATIONS FOR D END OF BREAKWATER. 
WEVE CHARACTERISTICS 
ANGLE 135-00 DF&Ges PemtOD 4e00 SE@es SPEEDZO 600 F/G. 
AT POINT X =25200 Y= 30.00 
DUeewecerIkRAVEL IS 180.17 DEGesGeae Or DIFFRACTION 1s. CereeG? 
WAVE CHARACTERISTICS 
PPYOLE 135.00 PFG.> PE D 4e00 SECes SPEEDO ZOROOP F | 7 aeG. 
AT POINT X =26.990 Y= 32.209 
DIRMOP TRAVEL IS 155.00 DEGes@eer OF DIFPRACTION 1S Os2466 
WAVE CHARACTERISTICS 
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ANGLE 135.200 DFGe» PE] OD 4eO00 SFCe>s APR ED2 0 «00 6F |S Ge 
AT POINT ..%4—~=2 7a GGe =. 323605 
DIR OF TRAVEL IS 141+20 DFG@ieeEF OF DIFERACTIONS IS O22 
WAVE CGHARAGTERISTICS 
SHEE 135.00 Pears PERIOD 4eGOGMoE Ges SPEEVZO<600. Fly 
eee OUNT X =2ZieS0meY= 3350 
DIR OF TRAVEL IS 13599 DEGuai@eeRe OF DIFFRAGRION SIS 3a 3c 
WAVE CHARACTERISTICS 
ANGLF 135200 DFGes PERIOD 4 V0 SSECes SPEEDZ@ COUR it aa. 
DIR OF TRAVEL- 1S 9132-59 DEGsesGeeepeece DIFFRACT DON MR SReles te >o 
DIFFRACTION CALCULATIONS FOR BREAKWATER GAPe 
BASED ON IMAGINARY GAP ASSUMED MORMAL TO WAVE RAYe 
IMAGINARY GAP FROM 2020003212000 TO 2220003212000 
WAVE CHARACTERISTICS 
ANGLE 90.90 DFGes PE) OD 4e00 SECes SeeeD20<«00 Fivsea. 
pepeeOINT X =271-90 Y= 26-00 
DIRT ORR AVE. TS 90-00 DEGesCOEF OF DlRB ee TION 1S) 20.c 473 
WAVE CHARACTERISTICS 
ANGLE 90.00 DEG. s PER LO Eel 0 SEG <:5 SPEEDZ0s:00> Fl 7sbG. 
AiePOINT X =23-00 Y= 25-409 
DIRSOR RAVE LTS 77¢82 D&GesCOFEF OF DIFFRAGTION 135 7Oe2266 
Wey eeeCrAeAGTER TST ICS 
ANGLF 90.00 DFGes PFRIOD 4000 SECes SPEED2Z0.600 Fi7 asc. 
Pee OINT X =20.50 Y= 22.09 
DikwObeIRAVEL IS 94629 DEGesCOEF OF DIFFRACTION TS re. see 
FLIMINATED BREAKWATERS 
BREAKWATER FROM 10.005 10-00 TO 20.00; 10200 
DIFFRACTION CALCULATIONS FOR A END OF BREAKWATER. 
WAVE CHARACTERISTICS 
ANGLE 2600 DEG. = PERLOG tive 0 Oo Gans SPEEDZO.00. F 17 Sec. 
AT POINT X = 9-90 Y= 15-00 
DER OF “TRAVEL 1S 84.98 DEG.sCOEF OF DIFFRAGTTONIS Osco” 
WAVE CHARACTERISTICS 
ANGLF oT .00 DRGs. PERTOD Lee SE Gees SPEEDZOIOO Jb ty ohGe 
Pome INT X =10.00 Y= 15.02 
DIR OF TRAVEL IS 84628 DEGes COEF OF DIFFRACTION Seco 
WAVE CHARACTERISTICS 
ANGLE 90.00 DFGes PERIOD 4ePO SECe>» SPEEDZOSO OCF l7 Sece 
DIR OF TRAVEL IS 83056 DEGesGQOEF CHEUDIFERRACTION 16° Cs4a%ce 
PigeeeeCTION CALCULATIONS FOR B END OF BREAKWATER. 
WAVE CHARACTERISTICS 
mG F SOs moO Ges PEKTOD Lg CD a Gears SPEEDZO0.0OCME I ste 
AT POINT X =192990 Y= 15.00 
DTR CFS@TRAVEL IS 96-44 DFG. sGORR oO, DIFFRACTION® IS O.4eo6e 
WAVE CHARACTERISTICS 
ANGLE 90.09 DFGes PERIOD 4) OW ae Ge = SPEEDZOsO0 F iva: Gs 
Pee OINT X =2)69@9 Y= 15.90 
DTineOgeITRAVEL 1S 950/22 DEGesGee, OF DIFFRACTION #S O68 5260 
WAVE CHARACTERISTICS 
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AN 


GRID 


OMNI AMF WNHY 2 


23 


OMNIA TDF NN HH 2 


SPEEDZ0 00056 17 Sea. 


085.569 


pear 
40264 
4182 
42610 
“lee 
Se 2 
36061 
34656 
3 ieee 
31283 
29643 
28294 
2S 
Cee 
260¢64 
24698 
Zaed 
2 Wen) 
20430 
18.57 
Loe 55 
ie ae gs. 
Ceo 
3010 


030 


DET 
4496 
43294 
42282 
42052 
42064 
4215 
40092 
Ses 
3409] 
Be 79 
52 a 
Silene 2 
2 eee 


Sie 90.0U DEGes Pee Ol) 4edO0 SECes 
At POINT X:=20 . eee 1 OD 
DiR OF Reve Is 950%? DE Gem@e@er OF Dipper h) Ohmis 
FLIMINATED RBREAKWATERS 
es 5 Oi 2B 40 WITH 300C29C FOOT SQUARES. 
ANALYTICAL DATA FROM X= TOTO 50 ANDY = 23) le 40 
TOPO DATA FROM X= Pe 50 AND Y= ea 0 32 
feet CASE 31 PERIOD = 4GeOSECe ANGLE= 302¢0DEGeX= 5-00Y=11219 
COme SHOAL 
X yi A G BETA RFTN COEF CURR 5 
5.0800 114100 30.000 20.42 1-CO9 Iii. 00 600 O0TAeA. Oe 
Oe(31 126102 30.159 20043 1-009 OFBReee.oe 0.0072 aeee. 
e459 135209 306291 20644 12-024 OsSeeGarIeC0CO 0460009 2500 
mee SS 146119 306392 20653 165445 Ceeeeeee “0.0017 2.00 
Me 208 156135 206527 20639 146076 W8ES ee. 0232.22.00 
Meo2O 1Ge15be40.9078 202.36 1lell2 06948 eet. 079 2.00 
ie oees) o l7el91 eee oO 20% 20 le eo 02°31 099 029033 2090 
Meme ee 1he7 ae 27146683 20626 12207? Geel? O06 Somee 46 2.00 
eet 42 1946292 326110 20621 1-426) 04691]. 0s9S 705 202 ee 
meee? 9 2C e267 226875 20609 14347 Dee62 0698  OsOlise 2.Co 
meee 4S 206640 34.6244 20696 1.4382 06351 04697 050729 50.50 
pez 06 2069S 446653 20604 124425 06838 0297 Oe0 leo 0.50 
fees! 2is 156 34.177 19499 1.476 Oe823 0697 O29771 0650 
22,092 2 ett 24.287] 19290 le Sa 02806 0297 O eke? 77 Oe 50 
ee O0 216767 3526786 19s74% 12618 OeveG 0696 Osb36B 0659 
Meme Os 226063 376001 19655 1el16 Oneesg 0295 O40487 Oeae 
fee 8 OP 26 210 286622 19624 124840 O86 F437 0495 020657 O6e0 
PeeooS 222689 402761 19e09 22-007? Os 707 0694 024.0840 9.50 
Oem 53 22.028 43.498 1875 2e216 O4672 0693 041087 0.450 
Bee 236380 470137 18225 22504 06632 0293 041493 0-50 
Meme 2361/62 5263279 Ive3! 26905 Ose OSS? 062710 Us5e 
meeO04 24,179 60e7 57 1504] 346489 06535 04693 0.635924 0.50 
Bee 169 244.650 B00545 92459 44397 OeG77 1209 1261307 0450 
ee > TOPPED 
mee CASE 32 PERION = 4-e9SECe ANGLFE= 30-0DEGeX= 6200Y= 9 
COEF’ SHOAL 
X 7 A C BETA RFTN COE re Geuky S 
6.000 94 390 3902N09 Peo 1 O00 1 « &O @ fe 00 Ne WOOF See. 
ela? 196301 302.051 200645 126300 160900 1-00 9.0095 2.00 
Peto? 1163073 206104 20044 12-00% 146C0P 1600 0120002 2.00 
Meee? 12 A707 A001 22 20644 12-000 1)ePOC 1600 0001 2.00 
Mee 2 1346 710 380612727 20644 12-009 1200 1400 969000 2-00 
eet 2 1462715 206140 20444 1.007 064999 1.00 GsGQ 07 2.09 
WemmeoU 15672) 240.264 2064) 1.008 OO88864 1.90 0.0014 2.09 
We OG 1G 3? 3046471 20027 1e¢018 OsFS] O699 0.20022 2.09 
Mme 26 175350 320.706 20.31 1.030 0.965 0299 0.0015 2.00 
Peo Or las 206857 20628 1-040 OCB 0.98 04-0015 2200 
232260 19.402 22605 2024 ete 2. Ona 0298 020029 CeO 
Mee 9 (0 206441 316422 20612 16069 OsGGa7 0.98 040022 2.00 
2Poe6/3 2146490 71.844 20608 1.09] 96.957 0-97 969961 2-00 
meee 28 276 C2) 326245 2OewS 1elll Oe 0.9/7 0690071 1.209 
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26251 





28.364 
2a 202 
Beek 
20.936 
Boe. 496 
20.860 
31-266 
312665 
B2— 054 
32431 
e267 | 
Be 122 
43.39% 
332514 


2aGG 7 
iaenOa 
Dec 
232654 
232936 
Phi? 
2h oes 
246815 
Melos 
2 Sugete > 
25-804 
Deen 9 
26.599 
270984 


RAY STOPPED 


Mine! CASE 33 


ODMDANA NDE YNNK 2 


X 

7 we)00 
Sie 2 
10.463 
2.619% 
ee 2 
ie «65° 
Py, 384 
feel l2 
204840 
Pe DO 3 
244282 
262900 
pee 112 
292420 
31.2121 
was 1% 
63.605! 
24.496 
B65 3 | 
352746 
Boe 15/ 
Picie 965 
fe4969 
272364 
376749 
38.2114 
386444 
38.703 


y 
7a 
Se Od 
96702 

Lee 0 
1208 
l?geel 2 
en es 
14292] 
Lee 2? 
16.94? 
1729464 
Perm 2 71 
20.028 
214064 
272116 
Pate 1 80 
2 yet kD 
24626? 
24.2814 
25.092 
2503276 
222062 
Ze? © | 
ZOeao 
266 206 
262928 
2714704 
27.731 


mer a TOPPED 
SHORELINE COORDINATES 


th 
2 
3 


22 e667 
ea 1 
23686 
Be 9 OO 
346478 
G7 OO 
262793 
270846 
39.858 
426276 
45-2644 
Pale! 
624640 
Cee? 


PERIOD = 


A 
306 900 
ae 37 
20.066 
20.2108 
aCe 1 25 
a eS 
2008176 
a eee 
20.404 
20.790 
20.954 
31e¢214 
Se 2 2 
Pale ee 
AZ? 
3222/40 
@™el1 i> 
Bere 25 
24.205 
342868 
352/20 
364844 
aha & 
41.052 
45.200 
Zeer 3 
Oa 


LOMO 16413750 e8 0697 OeG 4. 
1969] lel7GpmOmry>? 0.96 (Oe iO? 
19.86 1.207% SOs? 0.96 / ee eo...) 
19.84 1e¢236 0.899 0.96 Oey Ogee 
1Ge77? le2l? “Qleeeo 04.96 Oe07048 02D 
19662 12327 “@Qileese. 0695 OF id. oo 
19641 146391 Oars 9.95 O8025400050 
19.208 16474 Veees 0294 020622 Oe 0 
Pere 2 1058) Oem 06593 OUT Se Ges 
18626 lel18 Olea 9o3. O«lV O44 
1700? 16896 Osfeeeeo.o? 0.1501] > Oes. 
he6ce> 20149 02682 0.92 04.2796 ee 
14.604 26574 04623 04695 045250 02450 
6689 30347 (0634 ee 508 0650 
ie USE Ce ANGLE= 30-e¢0DEGeX= 7e00Y= 7 
COEF SHOAL 
G BETA S72 COEF GURY S 
20047 16906 1.00C 1-00 2.0003 4400 
20647 12001 126000 1-0C 0-2-0003 2-200 
20646 16002 06°99 12.00 OscC0G se. 
2g oo 1<¢004 0-998 lees) 020063 22000 
20044 146905 04997 1-00 -9.0000 2.00 
20-4? e006). 797 1.00 =9.0002 2.00 
20044 126006 02°97 LOO QGe-0C0] 2000 
20042 Le20° 02995 1e00 eon 2e0U 
20 wR Tee O92 0699 «202 2400 
20 623 Witte 0.9088 Us9O' Cove oe. ce 
20630 1469029 0-986 04698 0460012 2-00 
20026 Weee 06983 0498 OeCGZ0 2e00 
2062] 18046 04978 0-98 0.0024 2.00 
20613 1426060 94971 904697 0260031 2.00 
Z0e6RC lef 06963 0597 Us O0eer za 
19 eG? le Tle 06 33-05 96ry I. Os) 2.08 
19.92 lell® “06746 Ue Ceme ie fel. OU 
196@6 1el42 O69 45 50s ee 8 ls 
19.78 let/@ Os 922 02 3Ge 2. 14> Ice 
19,€9 1.207 Oso? Dire 2. 200 (20 
19 cms 1e236 OO CRO ee Poa 6S Ue 
19.47 16280 6864 0275 0.02734 9.20 
19-234 1¢33°9 02-864 0.94 Oe 465 Oe 
19.06 146419 06839 04692 0-0706 904-59 
186450 1.532 VetOeeenss 041130 0.50 
17052 12695 OCTGGe eee? 041835 0.50 
15661 16944 Ome. o2? 062198 0850 
12641 26346 CpGerpe. 9S 0.6438 Cos0 
16000 162048 
22000) 164500 
22000 17-4090 


Sirol 


27691 
2Ee/6 
Lbe2> 
2 >e 9 
Le 
23092 
22023 
20e2 
Lae 
1696 
142/70 
165.9 
749 

eres 


20 


DPTH 
46.95 
46019 
446/72 
42094 
42eC4 
41-96 
4263 
404699 
38600 
Sea 
Cae 
Pere ek 
31-78 
Boe 2 ) 
Zia? & 
2/7260 
Cowie 
26010 
25027 
2404? 
23249 
22072 
21-80 
2 Oe 
17-51 
14-20 
19-30 

Seioe 





° = a 


uF 








(IN 


TERIM CARDS 


3-200C 
4.000 
5 eO0C 
Gao 
62000 
1@ O00 
8 2000 
8 e000 
9.2000 
10.000 
11-2900 
11.2000 
25000 
132000 
142000 
Wa 000 


Site) 


492909 
41-4000 
422000 
434-000 
44,000 
442000 
452000 
462000 
472000 
48-000 
492900 
492000 
50-0009 


BemeoHORELINE SEGMENT 


ALL DONFe 


172¢90C 
17e¢464 
172¢667 
18 e000 
182318 
18 e600 
192000 
192000 
192444 
Doers 
202000 
202000 
202400 
Diet os 
21-000 
26 000 


28240 

282520 
284643 
282800 
292000 
292000 
292400 
296556 
279006! 
268 15 
Z0eNCO 
302000 
320¢4N90 
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